********************************************************
****	SPIR DFSA NUTRITION PAPER 					****
****	Tables construction							****
****	July 2024									****
********************************************************

/*

This do-file contains replication code for Tables 1-7, and Annex Tables A1, A2, 
B1, C1, C2 and D1 in the paper 'Including Scalable Nutrition Interventions in 
a Graduation Model Program: Experimental Evidence from Ethiopia'.

It uses three datasets, the data in which were collected as part of the Strengthen 
the PSNP4 Institutions and Resilience trial (2018-2021) and that are here limited 
to variables relevant to the analysis in the paper from the three rounds of 
data collection.

The datasets are as follows:
	- spir_household data.dta		household-level dataset from all survey rounds 
									used for the main analysis
	- spir_child level data.dta		child level (1 row per child) dataset from all 
									survey rounds used for the anthropometrics analysis
	- spir_midline data.dta			child level (1 row per child) dataset from the 
									midline survey round used for Table D1, replication
									of results from the Food Policy 2022 paper

If you have any questions, reach out to Heleene Tambet, h.tambet@cgiar.org.

*/

	clear
	clear mata
	clear matrix
	set maxvar 100000
	set more off

	// store date
	global date = "`c(current_date)'"

	// set directories
	global location "`c(pwd)'"
	global datain 	"${location}"	// change if analysis datasets moved to a subfolder/elsewhere than the location of the do-file
	global tables 	"${location}"	// change if you want outputted tables to be saved elsewhere than the location of the do-file
	
	
	// main analysis dataset
	use "${datain}/spir_household data.dta", clear
	
	
	
//==============================================================================//
// 1. MODEL SETUP //
{	
	
	// calculate sample weights for the poor sample
		svyset _n
		preserve
			keep if poultrycash_hh==1 
			svy: tabulate N_star_poultry N_star_cash if T1==1
			global pr_wght_NS_P = e(Obs)[2,1] / e(N)
			global pr_wght_NS_C = e(Obs)[1,2] / e(N)
			svy: tabulate N_poultry N_cash if T2==1
			global pr_wght_N_P = e(Obs)[2,1] / e(N)
			global pr_wght_N_C = e(Obs)[1,2] / e(N)
			di "$pr_wght_NS_P  $pr_wght_N_P  $pr_wght_NS_C  $pr_wght_N_C"
		restore
			
	// store model
		global model_poor 	N_star_poultry N_star_cash T3 N_poultry N_cash woredafixed*
		
	// store footnotes
		global footnote_hh		"Notes: Estimates from the DFSA SPIR endline survey sample. Standard errors (in parentheses) are clustered at the kebele level. All models control for woreda level fixed effects, child's age and gender, whether the mother has any education, and whether the day before the interview was a fasting day. For outcomes that were measured at baseline, the model also controls for the baseline value. Asterisks indicate significance at the 10, 5 and 1 percent level and are calculated with respect to the standard errors."
		global footnote_poor	"Sample of poor households is determined by ranking households within kebeles based on land and asset index constructed at baseline where 10 poorest households out of 18 in each kebele are classified as 'poor'."
		global footnote_aste	"Notes: Estimates from the DFSA SPIR endline survey sample. All estimates are calculated following the method of Katz, Kling and Liebman (2007) and present the effect size relative to the standard deviation of the control arm. Asterisks indicate significance at the 10, 5 and 1 percent level and are calculated with respect to the standard errors."
		global footnote_attrit	"Notes: Estimates from the DFSA SPIR endline survey sample. Standard errors (in parentheses) are clustered at the kebele level. Asterisks indicate significance at the 10, 5 and 1 percent level and are calculated with respect to the standard errors. Models control for woreda level fixed effects. This table excludes the 40 households from the 10 kebeles that were excluded from the endline survey due to insecurity. For baseline characteristics that had a missing value, we replaced these values with 0s and controlled for it with a respective binary variable."
	
	// store formatting parameters
		global hlines 	10100000000010001000000010001
		global sdec		3\3\3\3\3\3\3\3\3\3\3\3\3\3\3\3\3\3\3\3\3\3\3\3\0

}




//==============================================================================//
// 2. MAIN TABLES
{	

	//==============================================================================//
	**# Table 1. BASELINE BALANCE //
	{
		
	preserve
	keep if poultrycash_hh==1	// poor only sample
	
		local rvars log_tot_exp_ae hh_nearest_town_km hh_size_r1 pf_age_r1 pf_educ_r1 pm_age_r1 pm_educ_r1 wdds_r1 mddw_r1 ///
			ch_agem_r1 ch_male_r1 ch_stunted_r1 ch_wasted_r1 ch_weight_measure_r1 ch_height_measure_r1 total_food_r1 ch_mmf_r1 ch_mad_r1
		local nrvars : word count `rvars'
		
		// store matrix for estimates
		matrix Dstd = J(`nrvars',10,.) 
		cap drop __0000*
		
		// Ns
			tabstat `rvars', stat(count) save
			mat N1=r(StatTotal)' 
			
		// means by group
			// Enhanced Nutrition w/ grants
				tabstat `rvars' if T1==1, stat(mean sd) save
				mat D1=r(StatTotal)' 
			// Enhanced Nutrition w/o grants
				tabstat `rvars' if T3==1, stat(mean sd) save
				mat D2=r(StatTotal)' 
			// Core Nutrition
				tabstat `rvars' if T2==1, stat(mean sd) save
				mat D3=r(StatTotal)'
			// Control
				tabstat `rvars' if T4==1, stat(mean sd) save
				mat D4=r(StatTotal)'

		local i = 1
		
		// regressions for tests of differences
			foreach var of local rvars {
			
				regress `var' T1 if (T1==1 | T4==1), vce(cluster kebele_code)
					test T1
					matrix Dstd[`i',1] = r(p)
					matrix Dstd[`i',2] =.
			
				regress `var' T3 if (T3==1 | T4==1), vce(cluster kebele_code)
					test T3
					matrix Dstd[`i',3] = r(p)
					matrix Dstd[`i',4] =.
					
				regress `var' T2 if (T2==1 | T4==1), vce(cluster kebele_code)
					test T2
					matrix Dstd[`i',5] = r(p)
					matrix Dstd[`i',6] =.
					
				regress `var' T1 if (T1==1 | T3==1), vce(cluster kebele_code)
					test T1
					matrix Dstd[`i',7] = r(p)
					matrix Dstd[`i',8] =.
					
				regress `var' T1 if (T1==1 | T2==1), vce(cluster kebele_code)
					test T1
					matrix Dstd[`i',9] = r(p)
					matrix Dstd[`i',10] =.
					
				local ++i
				
			}
			mat list Dstd
					
		// create matrix for table
			mat N2 = J(`nrvars',1,.)
			mat D=N1,N2,D1,D2,D3,D4,Dstd
			mat list D

		// create table
		frmttable using "${tables}/all tables_${date}.doc", replace landscape statmat(D) varlabels substat(1) /// 
			sdec(0,3,3,3,3,3) ctitles("","N","{Mean and standard deviation}","","","","{p-value}","","","","" \ "","", ///
				"{EN with grants}""{EN without grants}","{CN with grants}","{Control}", ///
				"{EN + grants vs Control}","{EN, no grants vs Control}","{CN + grants vs Control}","{EN + grants vs EN, no grants}","{EN + grants vs CN + grants}") ///
			title("{Table 1: Balance in baseline characteristics}") multicol(1,3,4;1,7,5)  hlines(111{0}1) vlines(0 1 1 0 0 0 1 0 0 0 0) ///
			colw(50 4 4 4 4 4 4 4 4 4 4) titlfont(fs13 b) ctitlfont(fs10) rtitlfont(fs10) statfont(fs10) notefont(fs9) ///
			note(Notes: Estimates from the DFSA SPIR baseline survey sample. Standard deviations are in parentheses. P-value is from the test of difference of means between the treatment arms. Standard errors are clustered at the kebele level.)

	restore
	
	}
	
	
	//==============================================================================//
	**# Table 2: ANTHROPOMETRY: extremely poor //
	{
	
	preserve
	
	// anthropometrics data (1 row per child)
		use "$datain/spir_child level data.dta", clear
	
	// run the model & output the table
		keep if poultrycash_hh==1
		keep if !mi(ch_haz_r3) & !mi(ch_stunted_r3) & !mi(ch_whz_r3) & !mi(ch_wasted_r3) & !mi(ch_waz_r3) & !mi(ch_underwght_r3) & !mi(ch_muac_r3) & !mi(ch_sam_r3)
		
	// define matrices for statmat
		local rvars ch_haz ch_stunted ch_whz ch_wasted ch_waz ch_underwght ch_muac ch_sam	
		local i = 0
		local c : word count `rvars'
		local 2c = `c'*2
		clear matrix
		mat	N				= J(1,`2c',0)		// sample N
		mat	M				= J(1,`2c',0)		// control means
		  
	// matrices for model with extended controls
		mat C1_beta  		= J(5,`2c',0)		// coefficients for beta1 - beta9
		forval t = 1/4 {
			mat C1_lincom`t'= J(1,`2c',0)		// effects for linear combinations
		}
		forval t = 1/4 {
			mat C1_test`t'	= J(1,`2c',0)		// tests of equality
		}
		mat C1_pvalues		= J(11,`2c',0)		// p-values
		mat C1_stars  		= J(11,`2c',0)		// annotated coefficients with ***
							 
	// run loop to create matrices as individual inputs for statmat
		local z = 1
		
		forval i = 1/`c' {    
			
			tokenize `rvars'
			local rnames "`rvars'"
			 
			local end ``i''_r3
			
			// means (matrix M)
				mean `end' if `end'!=. & T4==1 & ch_agem_r3 >= 6 & ch_agem_r3 < 36
					mat M[1,`z'] = e(b)
					mat M[1,`z'+1] = .

			 // regression coefficients
				/*	Poultry X EN				β1
					Cash X EN					β2
					EN w/o grants			 	β3
					Poultry X CN				β4
					Cash x CN					β5		*/
				regress `end' $model_poor ch_agem_r3 ch_male_r3 pf_educ_r1 ///
						if ch_agem_r3 >= 6 & ch_agem_r3 < 36, vce(cluster kebele)
					mat P = r(table)
					mat N[1,`z'] = e(N)
					mat N[1,`z'+1] = .
					mat C1_beta[1,`z']=(_b[N_star_poultry],_se[N_star_poultry]\_b[N_star_cash],_se[N_star_cash]\ ///
						_b[T3],_se[T3]\_b[N_poultry],_se[N_poultry]\_b[N_cash],_se[N_cash])
			
			// p-values
				forval x = 1/5 {
					mat C1_pvalues[`x',`z'] = P["pvalue",`x']
				}
				
			// linear combination of betas
				/*	Linear combination of EN+grants		Weighted average: .5 * β1 + .5 * β2
					Linear combination of CN+grants		Weighted average: .5 * β4 + .5 * β5					*/					
				lincom N_star_poultry*$pr_wght_NS_P + N_star_cash*$pr_wght_NS_C
					mat C1_lincom1[1,`z'] = r(estimate), r(se)		// EN with grants
					mat C1_pvalues[6,`z'] = r(p)
				lincom N_poultry*$pr_wght_N_P + N_cash*$pr_wght_N_C
					mat C1_lincom2[1,`z'] = r(estimate), r(se)		// CN with grants
					mat C1_pvalues[7,`z'] = r(p)
	
			// test equality
				test N_star_poultry*$pr_wght_NS_P + N_star_cash*$pr_wght_NS_C = N_poultry*$pr_wght_N_P + N_cash*$pr_wght_N_C
					mat C1_test1[1,`z'] = r(p),.
				test N_star_poultry*$pr_wght_NS_P + N_star_cash*$pr_wght_NS_C = T3
					mat C1_test2[1,`z'] = r(p),.
				test N_star_poultry = N_star_cash
					mat C1_test3[1,`z'] = r(p),.
				test N_poultry = N_cash
					mat C1_test4[1,`z'] = r(p),.

			// p-values
				local y = 1
				forval x = 8/11 {
					mat C1_pvalues[`x',`z'] = C1_test`y'[1,`z']
					local y = `y' + 1
				}
								   
			 // annotate arm coefficients with sig. stars
				forval x = 1/11 {
					mat C1_stars[`x',`z'] = ((C1_pvalues[`x',`z']<=0.10) + (C1_pvalues[`x',`z']<=0.05) + (C1_pvalues[`x',`z']<=0.01)),0
				}
			mat drop P          // dropping matrix with regression results
			local z = `z'+2
			
		}         	
	
	// aggregate matrices into full statmat as input for frmttable         
		mat D = C1_beta\C1_lincom1\C1_lincom2\C1_test1\C1_test2\C1_test3\C1_test4\M\N		
		mat list D
		
	frmttable using "${tables}/all tables_$date.doc", addtable landscape statmat(D) substat(1) annotate(C1_stars) asymbol(*,**,***) ///
		hlines($hlines) sdec($sdec) ///
		title("Table 2: Anthropometrics, children 6-35 months of age") nobl ///
		basefont(fs10) titlfont(fs13 b) ctitlfont(fs10) rtitlfont(fs10) notefont(fs9) colwidth(28 10 10 10 10 10 10 10 10 10 10) ///
		rtitle("EN x Poultry"\""\"EN x Cash"\""\"EN without grants"\""\"CN x Poultry"\""\"CN x Cash"\""\ ///
			"Average effect of EN with grants"\""\"Average effect of CN with grants"\""\ ///
			"Test: EN with grants = CN with grants"\""\"Test: EN with grants = EN without grants"\""\"Test: EN x Poultry = EN x Cash"\""\"Test: CN x Poultry = CN x Cash"\""\"Mean of control"\""\"N")  ///
		ctitles("","(1)","(2)","(3)","(4)","(5)","(6)","(7)","(8)"\"","Height-for-age z-score (HAZ)","Proportion stunted (HAZ<-2SD)","Weight-for-height z-score (WHZ)", ///
			"Proportion wasted (WHZ<-2SD)","Weight-for-age z-score (WAZ)","Proportion underweight (WAZ<-2SD)","Mid-upper arm circumference (MUAC)", ///
			"Proportion of severe acute malnutrition (MUAC<11.5 cm)") ///
		note("Notes: Estimates from the DFSA SPIR endline survey for children of age 6 to 35 months from the subsample of extremely poor households. Standard errors (in parentheses) are clustered at the kebele level. All models control for woreda level fixed effects, child's sex and age in months, and whether child's mother has any education. Asterisks indicate significance at the 10, 5 and 1 percent level and are calculated with respect to the standard errors. ${footnote_poor}") ///
		coljust(l;c)	
			
	cap mat drop D N M I_0 C1_beta C1_lincom1 C1_lincom2 C1_lincom3 C1_lincom4 C1_test1 C1_test2 C1_test3 C1_test4 C1_test5 C1_test6 C1_test7 C1_pvalues C1_stars
	restore
	
	//============================================================================//		
	}

	
	//==============================================================================//
	**# Table 3: MOTHER'S ACCESS TO HEALTHCARE: extremely poor //
	{
		
	// replace missing BL outcomes with 0s + control for it
		foreach var in j7_hew_hda j7_hew_hda_home j7_fdemo_3mo j7_bcc_3mo ///
			wdds mddw wfoodgroup1 wfoodgroup2 wfoodgroup3 wfoodgroup4 wfoodgroup5 wfoodgroup6 wfoodgroup7 wfoodgroup8 wfoodgroup9 wfoodgroup10 ///
			knowledge_f knowledge_1 knowledge_3 knowledge_5 knowledge_11 knowledge_13 knowledge_19 knowledge_21 {
				gen miss_`var'_r1 = missing(`var'_r1)
				replace `var'_r1 = 0 if missing(`var'_r1) & !missing(`var'_r3)
		}
			
	preserve 
	
	keep if !mi(j7_hew_hda_r3) & !mi(j7_hew_hda_home_r3) & !mi(j7_fdemo_3mo_r3) & ///
		!mi(j7_bcc_3mo_r3) & !mi(j7_cltsh_r3) & !mi(j7_odf_r3) & !mi(m2_topic_2a_f_r3)
	keep if poultrycash_hh==1
	
	// define matrices for statmat
		local rvars j7_hew_hda j7_hew_hda_home j7_fdemo_3mo j7_bcc_3mo j7_cltsh j7_odf m2_topic_2a_f		
		local i = 0
		local c : word count `rvars'
		local 2c = `c'*2+2
		clear matrix
		mat	N				= J(1,`2c',0)		// sample N
		mat	M				= J(1,`2c',0)		// control means
		  
	// matrices for model with extended controls
		mat C1_beta  		= J(5,`2c',0)		// coefficients for beta1 - beta9
		forval t = 1/4 {
			mat C1_lincom`t'= J(1,`2c',0)		// effects for linear combinations
		}
		forval t = 1/4 {
			mat C1_test`t'	= J(1,`2c',0)		// tests of equality
		}
		mat C1_pvalues		= J(11,`2c',0)		// p-values
		mat C1_stars  		= J(11,`2c',0)		// annotated coefficients with ***
							 
	// run loop to create matrices as individual inputs for statmat
		local z = 3
		
		forval i = 1/`c' {    
			
			tokenize `rvars'
			local rnames "`rvars'"
			 
			local end ``i''_r3
			capture confirm variable ``i''_r1
			if !_rc {
				local base 	``i''_r1  miss_``i''_r1
			}
			else {
				local base 	""
			}
			
			// means (matrix M)
				mean `end' if `end'!=. & T4==1
					mat M[1,`z'] = e(b)
					mat M[1,`z'+1] = .

			 // regression coefficients
				regress `end' $model_poor pf_educ_r1 `base', vce(cluster kebele_code)
					mat P = r(table)
					mat N[1,`z'] = e(N)
					mat N[1,`z'+1] = .
					mat C1_beta[1,`z']=(_b[N_star_poultry],_se[N_star_poultry]\_b[N_star_cash],_se[N_star_cash]\ ///
						_b[T3],_se[T3]\_b[N_poultry],_se[N_poultry]\_b[N_cash],_se[N_cash])
			
			// p-values
				forval x = 1/5 {
					mat C1_pvalues[`x',`z'] = P["pvalue",`x']
				}
				
			// linear combination of betas				
				lincom N_star_poultry*$pr_wght_NS_P + N_star_cash*$pr_wght_NS_C
					mat C1_lincom1[1,`z'] = r(estimate), r(se)		// EN with grants
					mat C1_pvalues[6,`z'] = r(p)
				lincom N_poultry*$pr_wght_N_P + N_cash*$pr_wght_N_C
					mat C1_lincom2[1,`z'] = r(estimate), r(se)		// CN with grants
					mat C1_pvalues[7,`z'] = r(p)
					
			// test equality
				test N_star_poultry*$pr_wght_NS_P + N_star_cash*$pr_wght_NS_C = N_poultry*$pr_wght_N_P + N_cash*$pr_wght_N_C
					mat C1_test1[1,`z'] = r(p),.
				test N_star_poultry*$pr_wght_NS_P + N_star_cash*$pr_wght_NS_C = T3
					mat C1_test2[1,`z'] = r(p),.
				test N_star_poultry = N_star_cash
					mat C1_test3[1,`z'] = r(p),.
				test N_poultry = N_cash
					mat C1_test4[1,`z'] = r(p),.

			// p-values
				local y = 1
				forval x = 8/11 {
					mat C1_pvalues[`x',`z'] = C1_test`y'[1,`z']
					local y = `y' + 1
				}
								   
			 // annotate arm coefficients with sig. stars 
				forval x = 1/11 {
					mat C1_stars[`x',`z'] = ((C1_pvalues[`x',`z']<=0.10) + (C1_pvalues[`x',`z']<=0.05) + (C1_pvalues[`x',`z']<=0.01)),0
				}
			mat drop P          // dropping matrix with regression results
			local z = `z'+2
		  }  
		  
			// add ASTE 
				local z = 1
				local full_set j7_hew_hda_r3 j7_hew_hda_home_r3 j7_fdemo_3mo_r3 j7_bcc_3mo_r3 j7_cltsh_r3 j7_odf_r3 m2_topic_2a_f_r3
				avg_effect `full_set', x(N_star_poultry N_star_cash T3 N_poultry N_cash woredafixed* pf_educ_r1 `base') ///
					effect(N_star_poultry N_star_cash T3 N_poultry N_cash) control(T4==1) cluster(kebele_code) keepmissing
								
					// add coefficients + p-values to matrix cells
					local i = 1
					foreach var in N_star_poultry N_star_cash T3 N_poultry N_cash {
						
						// coefficient
							mat C1_beta[`i',`z'] = (_b[ae_`var'],_se[ae_`var'])
						// p-value
							mat C1_pvalues[`i',`z'] = 2*ttail(e(N)-1,abs(_b[ae_`var']/_se[ae_`var']))
						
						local i = `i'+1
					}
						
					// mean + N
						mat M[1,`z'] = .
						mat M[1,`z'+1] = .
						mat N[1,`z'] = e(N)
						mat N[1,`z'+1] = .
				
					// linear combination of betas				
						lincom _b[ae_N_star_poultry]*$pr_wght_NS_P + _b[ae_N_star_cash]*$pr_wght_NS_C
							mat C1_lincom1[1,`z'] = r(estimate), r(se)		// EN with grants
							mat C1_pvalues[6,`z'] = r(p)
						lincom _b[ae_N_poultry]*$pr_wght_N_P + _b[ae_N_cash]*$pr_wght_N_C
							mat C1_lincom2[1,`z'] = r(estimate), r(se)		// CN with grants
							mat C1_pvalues[7,`z'] = r(p)
		
					// test equality
						test _b[ae_N_star_poultry]*$pr_wght_NS_P + _b[ae_N_star_cash]*$pr_wght_NS_C = _b[ae_N_poultry]*$pr_wght_N_P + _b[ae_N_cash]*$pr_wght_N_C
							mat C1_test1[1,`z'] = r(p),.
							mat C1_pvalues[8,`z'] = C1_test1[1,`z']
						test _b[ae_N_star_poultry]*$pr_wght_NS_P + _b[ae_N_star_cash]*$pr_wght_NS_C = _b[ae_T3]
							mat C1_test2[1,`z'] = r(p),.
							mat C1_pvalues[9,`z'] = C1_test2[1,`z']
						test _b[ae_N_star_poultry] = _b[ae_N_star_cash]
							mat C1_test3[1,`z'] = r(p),.
							mat C1_pvalues[10,`z'] = C1_test3[1,`z']
						test _b[ae_N_poultry] = _b[ae_N_cash]
							mat C1_test4[1,`z'] = r(p),.
							mat C1_pvalues[11,`z'] = C1_test4[1,`z']
						
					 // annotate arm coefficients with sig. stars
						forval x = 1/11 {
							mat C1_stars[`x',`z'] = ((C1_pvalues[`x',`z']<=0.10) + (C1_pvalues[`x',`z']<=0.05) + (C1_pvalues[`x',`z']<=0.01)),0
						}	
	
	// aggregate matrices into full statmat as input for frmttable         
		mat D = C1_beta\C1_lincom1\C1_lincom2\C1_test1\C1_test2\C1_test3\C1_test4\M\N		
		mat list D
		
	frmttable using "${tables}/all tables_$date.doc", addtable landscape statmat(D) substat(1) annotate(C1_stars) asymbol(*,**,***) ///
		hlines($hlines) sdec($sdec) ///
		title("Table 3: Mother's access to healthcare") nobl ///
		basefont(fs10) titlfont(fs13 b) ctitlfont(fs10) rtitlfont(fs10) notefont(fs9) colwidth(28 10 10 10 10 10 10 10 10 10 10) ///
		rtitle("EN x Poultry"\""\"EN x Cash"\""\"EN without grants"\""\"CN x Poultry"\""\"CN x Cash"\""\ ///
			"Average effect of EN with grants"\""\"Average effect of CN with grants"\""\ ///
			"Test: EN with grants = CN with grants"\""\"Test: EN with grants = EN without grants"\""\"Test: EN x Poultry = EN x Cash"\""\"Test: CN x Poultry = CN x Cash"\""\ ///
			"Mean of control"\""\"N")  ///
		ctitles("","(1)","(2)","(3)","(4)","(5)","(6)","(7)","(8)"\"","Average Standardized Treatment Effect: Access to healthcare", ///
			"Has had contact with a HEW or a HDA in past 3 months","Has been visited by a HEW or a HDA at home in past 3 months", ///
			"Has attended a food demonstration in her community in last 3 months","Has attended a BCC session in past 3 months", ///
			"Has attended a community led total sanitation and hygiene (CLTSH) event","Has participated in a Open Defecation Free (ODF) event", ///
			"Has discussed nutrition topics at VESA meetings in past 12 months") ///
		note("Notes: Estimates from the DFSA SPIR endline survey for the subsample of extremely poor households. Standard errors (in parentheses) are clustered at the kebele level. All models control for woreda level fixed effects and whether mother has any education. Asterisks indicate significance at the 10, 5 and 1 percent level and are calculated with respect to the standard errors. ${footnote_poor}") ///
		coljust(l;c)	
			
	cap mat drop D N M I_0 C1_beta C1_lincom1 C1_lincom2 C1_lincom3 C1_lincom4 C1_test1 C1_test2 C1_test3 C1_test4 C1_test5 C1_test6 C1_test7 C1_pvalues C1_stars
	restore

	}
	
	
	//==============================================================================//
	**# Table 4. CHILD HEALTH HISTORY: extremely poor //
	{
	preserve 
	
	keep if eic_age_r3 >= 6 & eic_age_r3 < 24
	keep if !mi(ch_vita_r3) & !mi(ch_weight_measure_r3) & !mi(ch_height_measure_r3) & !mi(ch_muac_measure_r3) & !mi(ch_underweight_r3) & !mi(ch_cpnp_r3)
	keep if poultrycash_hh==1
	
	// define matrices for statmat
		local rvars ch_vita ch_weight_measure ch_height_measure ch_muac_measure ch_underweight ch_cpnp	
		local i = 0
		local c : word count `rvars'
		local 2c = `c'*2 + 2
		clear matrix
		mat	N				= J(1,`2c',0)		// sample N
		mat	M				= J(1,`2c',0)		// control means
		  
	// matrices for model with extended controls
		mat C1_beta  		= J(5,`2c',0)		// coefficients for beta1 - beta9
		forval t = 1/4 {
			mat C1_lincom`t'= J(1,`2c',0)		// effects for linear combinations
		}
		forval t = 1/4 {
			mat C1_test`t'	= J(1,`2c',0)		// tests of equality
		}
		mat C1_pvalues		= J(11,`2c',0)		// p-values
		mat C1_stars  		= J(11,`2c',0)		// annotated coefficients with ***
							 
	// run loop to create matrices as individual inputs for statmat
		local z = 3
		
		forval i = 1/`c' {    
			
			tokenize `rvars'
			local rnames "`rvars'"
			 
			local end ``i''_r3
			
			// means (matrix M)
				mean `end' if `end'!=. & T4==1
					mat M[1,`z'] = e(b)
					mat M[1,`z'+1] = .

			 // regression coefficients
				/*	Poultry X EN				β1
					Cash X EN					β2
					Poultry X CN			 	β3
					Cash X CN					β4
					T3/EN w/o grants			β5		*/
				regress `end' $model_poor eic_age_r3 eic_male_r3 pf_educ_r1, vce(cluster kebele_code)
					mat P = r(table)
					mat N[1,`z'] = e(N)
					mat N[1,`z'+1] = .
					mat C1_beta[1,`z']=(_b[N_star_poultry],_se[N_star_poultry]\_b[N_star_cash],_se[N_star_cash]\ ///
						_b[T3],_se[T3]\_b[N_poultry],_se[N_poultry]\_b[N_cash],_se[N_cash])
			
			// p-values
				forval x = 1/5 {
					mat C1_pvalues[`x',`z'] = P["pvalue",`x']
				}
				
			// linear combination of betas		
				lincom N_star_poultry*$pr_wght_NS_P + N_star_cash*$pr_wght_NS_C
					mat C1_lincom1[1,`z'] = r(estimate), r(se)		// EN with grants
					mat C1_pvalues[6,`z'] = r(p)
				lincom N_poultry*$pr_wght_N_P + N_cash*$pr_wght_N_C
					mat C1_lincom2[1,`z'] = r(estimate), r(se)		// CN with grants
					mat C1_pvalues[7,`z'] = r(p)
							
			// test equality
				test N_star_poultry*$pr_wght_NS_P + N_star_cash*$pr_wght_NS_C = N_poultry*$pr_wght_N_P + N_cash*$pr_wght_N_C
					mat C1_test1[1,`z'] = r(p),.
				test N_star_poultry*$pr_wght_NS_P + N_star_cash*$pr_wght_NS_C = T3
					mat C1_test2[1,`z'] = r(p),.
				test N_star_poultry = N_star_cash
					mat C1_test3[1,`z'] = r(p),.
				test N_poultry = N_cash
					mat C1_test4[1,`z'] = r(p),.

			// p-values
				local y = 1
				forval x = 8/11 {
					mat C1_pvalues[`x',`z'] = C1_test`y'[1,`z']
					local y = `y' + 1
				}
								   
			 // annotate arm coefficients with sig. stars
			 // pvalues are stored in columns 1,2 of matrix C1_pvalues 
				forval x = 1/11 {
					mat C1_stars[`x',`z'] = ((C1_pvalues[`x',`z']<=0.10) + (C1_pvalues[`x',`z']<=0.05) + (C1_pvalues[`x',`z']<=0.01)),0
				}
			mat drop P          // dropping matrix with regression results
			local z = `z'+2
		  }        
		  
			// add ASTE 
				local z = 1
				local full_set ch_vita_r3 ch_weight_measure_r3 ch_height_measure_r3 ch_muac_measure_r3 ch_underweight_r3 ch_cpnp_r3
				avg_effect `full_set', x(N_star_poultry N_star_cash T3 N_poultry N_cash woredafixed* eic_age_r3 eic_male_r3 pf_educ_r1) ///
					effect(N_star_poultry N_star_cash T3 N_poultry N_cash) control(T4==1) cluster(kebele_code) keepmissing
								
					// add coefficients + p-values to matrix cells
					local i = 1
					foreach var in N_star_poultry N_star_cash T3 N_poultry N_cash {
						
						// coefficient
							mat C1_beta[`i',`z'] = (_b[ae_`var'],_se[ae_`var'])
						// p-value
							mat C1_pvalues[`i',`z'] = 2*ttail(e(N)-1,abs(_b[ae_`var']/_se[ae_`var']))
						
						local i = `i'+1
					}
						
					// mean + N
						mat M[1,`z'] = .
						mat M[1,`z'+1] = .
						mat N[1,`z'] = e(N)
						mat N[1,`z'+1] = .
				
					// linear combination of betas				
						lincom _b[ae_N_star_poultry]*$pr_wght_NS_P + _b[ae_N_star_cash]*$pr_wght_NS_C
							mat C1_lincom1[1,`z'] = r(estimate), r(se)		// EN with grants
							mat C1_pvalues[6,`z'] = r(p)
						lincom _b[ae_N_poultry]*$pr_wght_N_P + _b[ae_N_cash]*$pr_wght_N_C
							mat C1_lincom2[1,`z'] = r(estimate), r(se)		// CN with grants
							mat C1_pvalues[7,`z'] = r(p)
			
					// test equality
						test _b[ae_N_star_poultry]*$pr_wght_NS_P + _b[ae_N_star_cash]*$pr_wght_NS_C = _b[ae_N_poultry]*$pr_wght_N_P + _b[ae_N_cash]*$pr_wght_N_C
							mat C1_test1[1,`z'] = r(p),.
							mat C1_pvalues[8,`z'] = C1_test1[1,`z']
						test _b[ae_N_star_poultry]*$pr_wght_NS_P + _b[ae_N_star_cash]*$pr_wght_NS_C = _b[ae_T3]
							mat C1_test2[1,`z'] = r(p),.
							mat C1_pvalues[9,`z'] = C1_test2[1,`z']
						test _b[ae_N_star_poultry] = _b[ae_N_star_cash]
							mat C1_test3[1,`z'] = r(p),.
							mat C1_pvalues[10,`z'] = C1_test3[1,`z']
						test _b[ae_N_poultry] = _b[ae_N_cash]
							mat C1_test4[1,`z'] = r(p),.
							mat C1_pvalues[11,`z'] = C1_test4[1,`z']
						
					 // annotate arm coefficients with sig. stars
						forval x = 1/11 {
							mat C1_stars[`x',`z'] = ((C1_pvalues[`x',`z']<=0.10) + (C1_pvalues[`x',`z']<=0.05) + (C1_pvalues[`x',`z']<=0.01)),0
						}	

	
	// aggregate matrices into full statmat as input for frmttable         
		mat D = C1_beta\C1_lincom1\C1_lincom2\C1_test1\C1_test2\C1_test3\C1_test4\M\N		
		mat list D
		
	frmttable using "${tables}/all tables_$date.doc", addtable landscape statmat(D) substat(1) annotate(C1_stars) asymbol(*,**,***) ///
		hlines($hlines) sdec($sdec) ///
		title("Table 4: Child's health history") nobl ///
		basefont(fs10) titlfont(fs13 b) ctitlfont(fs10) rtitlfont(fs10) notefont(fs9) colwidth(28 10 10 10 10 10 10 10 10 10 10) ///
		rtitle("EN x Poultry"\""\"EN x Cash"\""\"EN without grants"\""\"CN x Poultry"\""\"CN x Cash"\""\ ///
			"Average effect of EN with grants"\""\"Average effect of CN with grants"\""\ ///
			"Test: EN with grants = CN with grants"\""\"Test: EN with grants = EN without grants"\""\"Test: EN x Poultry = EN x Cash"\""\"Test: CN x Poultry = CN x Cash"\""\"Mean of control"\""\"N")  ///
		ctitles("","(1)","(2)","(3)","(4)","(5)","(6)","(7)"\"","Average Standardized Treatment Effect: Access to health and nutrition services", ///
			"Received dose of Vitamin A in past 6 months","Child's weight was measured in past 3 months","Child's height was measured in past 3 months", ///
			"Child's MUAC was measured in past 3 months","Child identified as severely or moderately underweight in past 6 months", ///
			"Child participated in a 2-week cooking demonstration and feeding session (CPNP)") ///
		note("Notes: Estimates from the DFSA SPIR endline survey for children aged 6-13 months from the extremely poor subsample. Standard errors (in parentheses) are clustered at the kebele level. All models control for woreda level fixed effects, child's gender and age in months, and whether mother has any education. Asterisks indicate significance at the 10, 5 and 1 percent level and are calculated with respect to the standard errors. ${footnote_poor}") ///
		coljust(l;c)	
			
	cap mat drop D N M I_0 C1_beta C1_lincom1 C1_lincom2 C1_lincom3 C1_lincom4 C1_test1 C1_test2 C1_test3 C1_test4 C1_test5 C1_test6 C1_test7 C1_pvalues C1_stars
	restore
	
	//============================================================================//		
	}
	
	
	//==============================================================================//
	**# Table 5: IYCF KNOWLEDGE OF WOMEN: extremely poor //
	{
	preserve
	
	keep if !mi(knowledge_f_r3) & !mi(knowledge_1_r3) & !mi(knowledge_3_r3) & !mi(knowledge_5_r3) & !mi(knowledge_11_r3) & !mi(knowledge_13_r3) & !mi(knowledge_19_r3) & !mi(knowledge_21_r3)
	keep if poultrycash_hh==1
	
	// define matrices for statmat
		local rvars knowledge_f knowledge_1 knowledge_3 knowledge_5 knowledge_11 knowledge_13 knowledge_19 knowledge_21
		local i = 0
		local c : word count `rvars'
		local 2c = `c'*2
		clear matrix
		mat	N				= J(1,`2c',0)		// sample N
		mat	M				= J(1,`2c',0)		// control means
		  
	// matrices for model with extended controls
		mat C1_beta  		= J(5,`2c',0)		// coefficients for beta1 - beta9
		forval t = 1/4 {
			mat C1_lincom`t'= J(1,`2c',0)		// effects for linear combinations
		}
		forval t = 1/4 {
			mat C1_test`t'	= J(1,`2c',0)		// tests of equality
		}
		mat C1_pvalues		= J(11,`2c',0)		// p-values
		mat C1_stars  		= J(11,`2c',0)		// annotated coefficients with ***
							 
	// run loop to create matrices as individual inputs for statmat
		local z = 1
		
		forval i = 1/`c' {    
			
			tokenize `rvars'
			local rnames "`rvars'"
			 
			local end ``i''_r3
			capture confirm variable ``i''_r1
			if !_rc {
				local base 	``i''_r1  miss_``i''_r1
			}
			else {
				local base 	""
			}
			
			// means (matrix M)
				mean `end' if `end'!=. & T4==1
					mat M[1,`z'] = e(b)
					mat M[1,`z'+1] = .

			 // regression coefficients
				/*	Poultry X EN				β1
					Cash X EN					β2
					Poultry X CN			 	β3
					Cash X CN					β4
					T3/EN w/o grants			β5		*/
				regress `end' $model_poor pf_educ_r1 `base', vce(cluster kebele_code)
					mat P = r(table)
					mat N[1,`z'] = e(N)
					mat N[1,`z'+1] = .
					mat C1_beta[1,`z']=(_b[N_star_poultry],_se[N_star_poultry]\_b[N_star_cash],_se[N_star_cash]\ ///
						_b[T3],_se[T3]\_b[N_poultry],_se[N_poultry]\_b[N_cash],_se[N_cash])
			
			// p-values
				forval x = 1/5 {
					mat C1_pvalues[`x',`z'] = P["pvalue",`x']
				}
				
			// linear combination of betas					
				lincom N_star_poultry*$pr_wght_NS_P + N_star_cash*$pr_wght_NS_C
					mat C1_lincom1[1,`z'] = r(estimate), r(se)		// EN with grants
					mat C1_pvalues[6,`z'] = r(p)
				lincom N_poultry*$pr_wght_N_P + N_cash*$pr_wght_N_C
					mat C1_lincom2[1,`z'] = r(estimate), r(se)		// CN with grants
					mat C1_pvalues[7,`z'] = r(p)
					
			// test equality
				test N_star_poultry*$pr_wght_NS_P + N_star_cash*$pr_wght_NS_C = N_poultry*$pr_wght_N_P + N_cash*$pr_wght_N_C
					mat C1_test1[1,`z'] = r(p),.
				test N_star_poultry*$pr_wght_NS_P + N_star_cash*$pr_wght_NS_C = T3
					mat C1_test2[1,`z'] = r(p),.
				test N_star_poultry = N_star_cash
					mat C1_test3[1,`z'] = r(p),.
				test N_poultry = N_cash
					mat C1_test4[1,`z'] = r(p),.

			// p-values
				local y = 1
				forval x = 8/11 {
					mat C1_pvalues[`x',`z'] = C1_test`y'[1,`z']
					local y = `y' + 1
				}
								   
			 // annotate arm coefficients with sig. stars
			 // pvalues are stored in columns 1,2 of matrix C1_pvalues 
				forval x = 1/11 {
					mat C1_stars[`x',`z'] = ((C1_pvalues[`x',`z']<=0.10) + (C1_pvalues[`x',`z']<=0.05) + (C1_pvalues[`x',`z']<=0.01)),0
				}
			mat drop P          // dropping matrix with regression results
			local z = `z'+2
		  }         	
	
	// aggregate matrices into full statmat as input for frmttable         
		mat D = C1_beta\C1_lincom1\C1_lincom2\C1_test1\C1_test2\C1_test3\C1_test4\M\N		
		mat list D
		
	frmttable using "${tables}/all tables_$date.doc", addtable landscape statmat(D) substat(1) annotate(C1_stars) asymbol(*,**,***) ///
		hlines($hlines) sdec($sdec) ///
		title("Table 5: Mother's IYCF knowledge") nobl ///
		basefont(fs10) titlfont(fs13 b) ctitlfont(fs10) rtitlfont(fs10) notefont(fs9) colwidth(28 10 10 10 10 10 10 10 10 10 10) ///
		rtitle("EN x Poultry"\""\"EN x Cash"\""\"EN without grants"\""\"CN x Poultry"\""\"CN x Cash"\""\ ///
			"Average effect of EN with grants"\""\"Average effect of CN with grants"\""\ ///
			"Test: EN with grants = CN with grants"\""\"Test: EN with grants = EN without grants"\""\"Test: EN x Poultry = EN x Cash"\""\"Test: CN x Poultry = CN x Cash"\""\"Mean of control"\""\"N")  ///
		ctitles("","(1)","(2)","(3)","(4)","(5)","(6)","(7)","(8)"\"","Female IYCF knowledge score (0-7)", ///
			"Knows how long after birth should a baby start breastfeeding","Knows until what age a baby should be exclusively breastfed", ///
			"Knows what to do if a mother thinks her baby is not getting enough breast milk","Knows which foods are rich in vitamin A", ///
			"Knows the common problem with gruels given as first foods to babies","Knows how often a baby 6-23 months old should eat animal source foods", ///
			"Knows how often a child should be fed when sick") ///
		note("Notes: Estimates from the DFSA SPIR endline survey for the subsample of extremely poor households. Standard errors (in parentheses) are clustered at the kebele level. All models control for woreda level fixed effects and whether mother has any education. Asterisks indicate significance at the 10, 5 and 1 percent level and are calculated with respect to the standard errors. ${footnote_poor}") ///
		coljust(l;c)	
			
	cap mat drop D N M I_0 C1_beta C1_lincom1 C1_lincom2 C1_lincom3 C1_lincom4 C1_test1 C1_test2 C1_test3 C1_test4 C1_test5 C1_test6 C1_test7 C1_pvalues C1_stars
	restore
	
	//============================================================================//		
	}
	
	
	//==============================================================================//
	**# Table 6: DIETARY DIVERSITY FOR WOMEN: extremely poor //
	{
	preserve 
	
	keep if !mi(wdds_r3) & !mi(mddw_r3) & !mi(wfoodgroup1_r3) & !mi(wfoodgroup2_r3) & !mi(wfoodgroup3_r3) & !mi(wfoodgroup4_r3) & ///
		!mi(wfoodgroup5_r3) & !mi(wfoodgroup6_r3) & !mi(wfoodgroup7_r3) & !mi(wfoodgroup8_r3) & !mi(wfoodgroup9_r3) & !mi(wfoodgroup10_r3)
	keep if poultrycash_hh==1
	
	// define matrices for statmat
		local rvars wdds mddw wfoodgroup1 wfoodgroup2 wfoodgroup3 wfoodgroup4 wfoodgroup5 wfoodgroup6 wfoodgroup7 wfoodgroup8 wfoodgroup9 wfoodgroup10 
		local i = 0
		local c : word count `rvars'
		local 2c = `c'*2
		clear matrix
		mat	N				= J(1,`2c',0)		// sample N
		mat	M				= J(1,`2c',0)		// control means
		  
	// matrices for model with extended controls
		mat C1_beta  		= J(5,`2c',0)		// coefficients for beta1 - beta9
		forval t = 1/4 {
			mat C1_lincom`t'= J(1,`2c',0)		// effects for linear combinations
		}
		forval t = 1/4 {
			mat C1_test`t'	= J(1,`2c',0)		// tests of equality
		}
		mat C1_pvalues		= J(11,`2c',0)		// p-values
		mat C1_stars  		= J(11,`2c',0)		// annotated coefficients with ***
							 
	// run loop to create matrices as individual inputs for statmat
		local z = 1
		
		forval i = 1/`c' {    
			
			tokenize `rvars'
			local rnames "`rvars'"
			 
			local end ``i''_r3
			capture confirm variable ``i''_r1
			if !_rc {
				local base 	``i''_r1  miss_``i''_r1
			}
			else {
				local base 	""
			}
			
			// means (matrix M)
				mean `end' if `end'!=. & T4==1
					mat M[1,`z'] = e(b)
					mat M[1,`z'+1] = .

			 // regression coefficients
				regress `end' $model_poor post_fast_r3 pf_educ_r1 `base', vce(cluster kebele_code)
					mat P = r(table)
					mat N[1,`z'] = e(N)
					mat N[1,`z'+1] = .
					mat C1_beta[1,`z']=(_b[N_star_poultry],_se[N_star_poultry]\_b[N_star_cash],_se[N_star_cash]\ ///
						_b[T3],_se[T3]\_b[N_poultry],_se[N_poultry]\_b[N_cash],_se[N_cash])
			
			// p-values
				forval x = 1/5 {
					mat C1_pvalues[`x',`z'] = P["pvalue",`x']
				}
				
			// linear combination of betas					
				lincom N_star_poultry*$pr_wght_NS_P + N_star_cash*$pr_wght_NS_C
					mat C1_lincom1[1,`z'] = r(estimate), r(se)		// EN with grants
					mat C1_pvalues[6,`z'] = r(p)
				lincom N_poultry*$pr_wght_N_P + N_cash*$pr_wght_N_C
					mat C1_lincom2[1,`z'] = r(estimate), r(se)		// CN with grants
					mat C1_pvalues[7,`z'] = r(p)
					
			// test equality
				test N_star_poultry*$pr_wght_NS_P + N_star_cash*$pr_wght_NS_C = N_poultry*$pr_wght_N_P + N_cash*$pr_wght_N_C
					mat C1_test1[1,`z'] = r(p),.
				test N_star_poultry*$pr_wght_NS_P + N_star_cash*$pr_wght_NS_C = T3
					mat C1_test2[1,`z'] = r(p),.
				test N_star_poultry = N_star_cash
					mat C1_test3[1,`z'] = r(p),.
				test N_poultry = N_cash
					mat C1_test4[1,`z'] = r(p),.

			// p-values
				local y = 1
				forval x = 8/11 {
					mat C1_pvalues[`x',`z'] = C1_test`y'[1,`z']
					local y = `y' + 1
				}
								   
			 // annotate arm coefficients with sig. stars
			 // pvalues are stored in columns 1,2 of matrix C1_pvalues 
				forval x = 1/11 {
					mat C1_stars[`x',`z'] = ((C1_pvalues[`x',`z']<=0.10) + (C1_pvalues[`x',`z']<=0.05) + (C1_pvalues[`x',`z']<=0.01)),0
				}
			mat drop P          // dropping matrix with regression results
			local z = `z'+2
		  }         	
	
	// aggregate matrices into full statmat as input for frmttable         
		mat D = C1_beta\C1_lincom1\C1_lincom2\C1_test1\C1_test2\C1_test3\C1_test4\M\N		
		mat list D
		
	frmttable using "${tables}/all tables_$date.doc", addtable landscape statmat(D) substat(1) annotate(C1_stars) asymbol(*,**,***) ///
		hlines($hlines) sdec($sdec) ///
		title("Table 6: Mother's dietary diversity") nobl ///
		basefont(fs10) titlfont(fs13 b) ctitlfont(fs10) rtitlfont(fs10) notefont(fs9) colwidth(28 10 10 10 10 10 10 10 10 10 10) ///
		rtitle("EN x Poultry"\""\"EN x Cash"\""\"EN without grants"\""\"CN x Poultry"\""\"CN x Cash"\""\ ///
			"Average effect of EN with grants"\""\"Average effect of CN with grants"\""\ ///
			"Test: EN with grants = CN with grants"\""\"Test: EN with grants = EN without grants"\""\"Test: EN x Poultry = EN x Cash"\""\"Test: CN x Poultry = CN x Cash"\""\"Mean of control"\""\"N")  ///
		ctitles("","(1)","(2)","(3)","(4)","(5)","(6)","(7)","(8)","(9)","(10)","(11)","(12)"\"", ///
			"Women's Dietary Diversity Score, WDDS (1-10)","Meets Minimum Dietary Diversity for Women, MDD-W" ///
			"Mother consumed tubers and grains","Mother consumed pulses", ///
			"Mother consumed nuts and seeds","Mother consumed dairy","Mother consumed meat, fish, poultry", ///
			"Mother consumed eggs","Mother consumed green leafy vegetables", ///
			"Mother consumed vitamin A rich fruits and vegetables","Mother consumed other vegetables", ///
			"Mother consumed other fruits")  ///
		note("Notes: Estimates from the DFSA SPIR endline survey for the subsample of extremely poor households. Standard errors (in parentheses) are clustered at the kebele level. All models control for woreda level fixed effects, whether mother has any education, and whether the day before the survey was a fast day. For outcomes that were measured at baseline, the model also controls for the baseline value. Asterisks indicate significance at the 10, 5 and 1 percent level and are calculated with respect to the standard errors. ${footnote_poor}") ///
		coljust(l;c)	
			
	cap mat drop D N M I_0 C1_beta C1_lincom1 C1_lincom2 C1_lincom3 C1_lincom4 C1_test1 C1_test2 C1_test3 C1_test4 C1_test5 C1_test6 C1_test7 C1_pvalues C1_stars
	restore
	
	//============================================================================//		
	}
	
	
	//==============================================================================//
	**# Table 7: DIETARY DIVERSITY FOR CHILDREN: extremely poor //
	{
	preserve
	
	keep if !mi(total_food_r3) & !mi(ch_mmf_r3) & !mi(foodgroup0_r3) & !mi(foodgroup1_r3) & !mi(foodgroup2_r3) & ///
		!mi(foodgroup3_r3) & !mi(foodgroup4_r3) & !mi(foodgroup5_r3) & !mi(foodgroup6_r3) & !mi(foodgroup7_r3)
	keep if poultrycash_hh==1
		
	// define matrices for statmat
		local rvars total_food ch_mdd ch_mmf foodgroup0 foodgroup1 foodgroup2 foodgroup3 foodgroup4 foodgroup5 foodgroup6 foodgroup7
		local i = 0
		local c : word count `rvars'
		local 2c = `c'*2
		clear matrix
		mat	N				= J(1,`2c',0)		// sample N
		mat	M				= J(1,`2c',0)		// control means
		  
	// matrices for model with extended controls
		mat C1_beta  		= J(5,`2c',0)		// coefficients for beta1 - beta9
		forval t = 1/4 {
			mat C1_lincom`t'= J(1,`2c',0)		// effects for linear combinations
		}
		forval t = 1/4 {
			mat C1_test`t'	= J(1,`2c',0)		// tests of equality
		}
		mat C1_pvalues		= J(11,`2c',0)		// p-values
		mat C1_stars  		= J(11,`2c',0)		// annotated coefficients with ***
							 
	// run loop to create matrices as individual inputs for statmat
		local z = 1
		
		forval i = 1/`c' {    
			
			tokenize `rvars'
			local rnames "`rvars'"
			 
			local end ``i''_r3
			
			// means (matrix M)
				mean `end' if `end'!=. & T4==1 & eic_age_r3 >= 6 & eic_age_r3 < 24
					mat M[1,`z'] = e(b)
					mat M[1,`z'+1] = .

			 // regression coefficients
				regress `end' $model_poor eic_age_r3 eic_male_r3 post_fast_r3 pf_educ_r1 ///
					if eic_age_r3 >= 6 & eic_age_r3 < 24, vce(cluster kebele_code)
					mat P = r(table)
					mat N[1,`z'] = e(N)
					mat N[1,`z'+1] = .
					mat C1_beta[1,`z']=(_b[N_star_poultry],_se[N_star_poultry]\_b[N_star_cash],_se[N_star_cash]\ ///
						_b[T3],_se[T3]\_b[N_poultry],_se[N_poultry]\_b[N_cash],_se[N_cash])
			
			// p-values
				forval x = 1/5 {
					mat C1_pvalues[`x',`z'] = P["pvalue",`x']
				}
				
			// linear combination of betas				
				lincom N_star_poultry*$pr_wght_NS_P + N_star_cash*$pr_wght_NS_C
					mat C1_lincom1[1,`z'] = r(estimate), r(se)		// EN with grants
					mat C1_pvalues[6,`z'] = r(p)
				lincom N_poultry*$pr_wght_N_P + N_cash*$pr_wght_N_C
					mat C1_lincom2[1,`z'] = r(estimate), r(se)		// CN with grants
					mat C1_pvalues[7,`z'] = r(p)
					
			// test equality
				test N_star_poultry*$pr_wght_NS_P + N_star_cash*$pr_wght_NS_C = N_poultry*$pr_wght_N_P + N_cash*$pr_wght_N_C
					mat C1_test1[1,`z'] = r(p),.
				test N_star_poultry*$pr_wght_NS_P + N_star_cash*$pr_wght_NS_C = T3
					mat C1_test2[1,`z'] = r(p),.
				test N_star_poultry = N_star_cash
					mat C1_test3[1,`z'] = r(p),.
				test N_poultry = N_cash
					mat C1_test4[1,`z'] = r(p),.

			// p-values
				local y = 1
				forval x = 8/11 {
					mat C1_pvalues[`x',`z'] = C1_test`y'[1,`z']
					local y = `y' + 1
				}
								   
			 // annotate arm coefficients with sig. stars
			 // pvalues are stored in columns 1,2 of matrix C1_pvalues 
				forval x = 1/11 {
					mat C1_stars[`x',`z'] = ((C1_pvalues[`x',`z']<=0.10) + (C1_pvalues[`x',`z']<=0.05) + (C1_pvalues[`x',`z']<=0.01)),0
				}
			mat drop P          // dropping matrix with regression results
			local z = `z'+2
		  }         	
	
	// aggregate matrices into full statmat as input for frmttable         
		mat D = C1_beta\C1_lincom1\C1_lincom2\C1_test1\C1_test2\C1_test3\C1_test4\M\N		
		mat list D
		
	frmttable using "${tables}/all tables_$date.doc", addtable landscape statmat(D) substat(1) annotate(C1_stars) asymbol(*,**,***) ///
		hlines($hlines) sdec($sdec) ///
		title("Table 7: Child's dietary diversity, children 6-23 months of age") nobl ///
		basefont(fs10) titlfont(fs13 b) ctitlfont(fs10) rtitlfont(fs10) notefont(fs9) colwidth(28 10 10 10 10 10 10 10 10 10 10 10) ///
		rtitle("EN x Poultry"\""\"EN x Cash"\""\"EN without grants"\""\"CN x Poultry"\""\"CN x Cash"\""\ ///
			"Average effect of EN with grants"\""\"Average effect of CN with grants"\""\ ///
			"Test: EN with grants = CN with grants"\""\"Test: EN with grants = EN without grants"\""\"Test: EN x Poultry = EN x Cash"\""\"Test: CN x Poultry = CN x Cash"\""\"Mean of control"\""\"N")  ///
		ctitles("","(1)","(2)","(3)","(4)","(5)","(6)","(7)","(8)","(9)","(10)","(11)"\"","Number of food groups (of 8) consumed", ///
			"Child meets meet minimum dietary diversity","Child meets minimum meal frequency", ///
			"Child consumed breast milk","Child consumed grains, roots or tubers","Child consumed legumes or nuts","Child consumed dairy","Child consumed fish or meat", ///
			"Child consumed eggs","Child consumed vitamin A rich fruits or vegetables","Child consumed other fruits or vegetables") ///
		note("Notes: Estimates from the DFSA SPIR endline survey for children aged 6-23 months from the subsample of extremely poor households. Standard errors (in parentheses) are clustered at the kebele level. All models control for woreda level fixed effects, child's gender and age in months, whether mother has any education, and whether the day before the survey was a fast day. Asterisks indicate significance at the 10, 5 and 1 percent level and are calculated with respect to the standard errors. ${footnote_poor}") ///
		coljust(l;c)	
			
	cap mat drop D N M I_0 C1_beta C1_lincom1 C1_lincom2 C1_lincom3 C1_lincom4 C1_test1 C1_test2 C1_test3 C1_test4 C1_test5 C1_test6 C1_test7 C1_pvalues C1_stars
	restore
	
	//============================================================================//		
	}

}


//==============================================================================//
// 3. ANNEX TABLES 
{
	
	//==============================================================================//
	**# A1: ATTRITION BY TREATMENT ARM //
	{		
	preserve 
	drop if inlist(kebele,"Selamgae","Shimila albeta","Bad","Cherey")	// drop the 4 kebeles that weren't surveyed due to unrest (see the paper)
	keep if poultrycash_hh==1 
	
	// missing BL controls - replace with 0 and control for it  
		foreach var in hh_size_r1 total_ch5_r1 log_tot_exp_ae hh_nearest_town_km ///
			pf_marital_1_r1 pf_age_r1 pf_educ_r1 pm_educ_r1 {
				gen m_bl_`var' = mi(`var')
				replace `var' = 0 if mi(`var')
			}
					
	// HOUSEHOLD LEVEL ATTRITION //
	
	// define matrices for statmat
		local rvars atr_hh_bl_el
		local i = 0
		local c : word count `rvars'
		local 3c = 4
		clear matrix
		mat	N				= J(1,`3c',0)		// sample N
		mat	M				= J(1,`3c',0)		// control means
		  
	// matrices for model with extended controls	
		mat C1_beta  		= J(10,`3c',0)		// coefficients for beta1 - beta3
		forval t = 1/3 {
			mat C1_test`t'	= J(1,`3c',0)		// tests of equality
		}
		mat C1_pvalues		= J(13,`3c',.)		// p-values
		mat C1_stars1  		= J(13,`3c',0)		// annotated coefficients with ***
					 
		local controls hh_size_r1 total_ch5_r1 log_tot_exp_ae hh_nearest_town_km ///
			pf_marital_1_r1 pf_age_r1 pf_educ_r1
			
	// run loop to create matrices as individual inputs for statmat
		local z = 1
		local x = 1 
			
			// means (matrix M)
				mean atr_hh_bl_el if atr_hh_bl_el!=. & T4==1
					mat M[1,`z'] = e(b)
					mat M[1,`z'+1] = .

			 // regression coefficients
				regress atr_hh_bl_el T1 T3 T2 woredafixed*, vce(cluster kebele_code)
					mat P = r(table)
					mat N[1,`z'] = e(N)
					mat N[1,`z'+1] = .
					mat C1_beta[1,`z']=(_b[T1],_se[T1]\_b[T3],_se[T3]\_b[T2],_se[T2]\ ///
						.,.\.,.\.,.\.,.\.,.\.,.\.,.)
	
			// p-values
				forval x = 1/3 {
					mat C1_pvalues[`x',`z'] = P["pvalue",`x']
				}
				
			// test equality
				test T1 = T2
					mat C1_test1[1,`z'] = r(p),.
				test T2 = T3
					mat C1_test2[1,`z'] = r(p),.
				test T1 = T3
					mat C1_test3[1,`z'] = r(p),.

			// p-values
				local y = 1
				forval x = 11/13 {
					mat C1_pvalues[`x',`z'] = C1_test`y'[1,`z']
					local y = `y' + 1
				}
								   
			 // annotate arm coefficients with sig. stars
				forval x = 1/13 {
					mat C1_stars1[`x',`z'] = ((C1_pvalues[`x',`z']<=0.10) + (C1_pvalues[`x',`z']<=0.05) + (C1_pvalues[`x',`z']<=0.01)),0
				}
				
			mat drop P          // dropping matrix with regression results
	
	// now with baseline controls added
	
	// run loop to create matrices as individual inputs for statmat
		local z = 3
		local x = 2 
			
			// means (matrix M)
				mean atr_hh_bl_el if atr_hh_bl_el!=. & T4==1
					mat M[1,`z'] = e(b)
					mat M[1,`z'+1] = .

			 // regression coefficients
				regress atr_hh_bl_el T1 T3 T2 `controls' m_bl_* woredafixed*, vce(cluster kebele_code)
					mat P = r(table)
					mat N[1,`z'] = e(N)
					mat N[1,`z'+1] = .
					mat C1_beta[1,`z']=(_b[T1],_se[T1]\_b[T3],_se[T3]\_b[T2],_se[T2]\ ///
						_b[hh_size_r1],_se[hh_size_r1]\_b[total_ch5_r1],_se[total_ch5_r1]\_b[log_tot_exp_ae],_se[log_tot_exp_ae]\ ///
						_b[hh_nearest_town_km],_se[hh_nearest_town_km]\_b[pf_marital_1_r1],_se[pf_marital_1_r1]\ ///
						_b[pf_age_r1],_se[pf_age_r1]\_b[pf_educ_r1],_se[pf_educ_r1])
	
			// p-values
				forval x = 1/10 {
					mat C1_pvalues[`x',`z'] = P["pvalue",`x']
				}
				
			// test equality
				test T1 = T2
					mat C1_test1[1,`z'] = r(p),.
				test T2 = T3
					mat C1_test2[1,`z'] = r(p),.
				test T1 = T3
					mat C1_test3[1,`z'] = r(p),.

			// p-values
				local y = 1
				forval x = 11/13 {
					mat C1_pvalues[`x',`z'] = C1_test`y'[1,`z']
					local y = `y' + 1
				}
								   
			 // annotate arm coefficients with sig. stars
				forval x = 1/13 {
					mat C1_stars1[`x',`z'] = ((C1_pvalues[`x',`z']<=0.10) + (C1_pvalues[`x',`z']<=0.05) + (C1_pvalues[`x',`z']<=0.01)),0
				}
				
			mat drop P          // dropping matrix with regression results
			
	// aggregate matrices into full statmat as input for frmttable         
		mat D1 = C1_beta\C1_test1\C1_test2\C1_test3\M\N	
	
	// CHILD LEVEL //
	
	// define matrices for statmat
		local rvars atr_ch_bl_el
		local i = 0
		local c : word count `rvars'
		local 3c = 4
		mat	N				= J(1,`3c',0)		// sample N
		mat	M				= J(1,`3c',0)		// control means
		  
	// matrices for model with extended controls	
		mat C1_beta  		= J(10,`3c',0)		// coefficients for beta1 - beta3
		forval t = 1/3 {
			mat C1_test`t'	= J(1,`3c',0)		// tests of equality
		}
		mat C1_pvalues		= J(13,`3c',.)		// p-values
		mat C1_stars2  		= J(13,`3c',0)		// annotated coefficients with ***
					 
		local controls hh_size_r1 total_ch5_r1 log_tot_exp_ae hh_nearest_town_km ///
			pf_marital_1_r1 pf_age_r1 pf_educ_r1
			
	// run loop to create matrices as individual inputs for statmat
		local z = 1
		local x = 1 
			
			// means (matrix M)
				mean atr_ch_bl_el if atr_ch_bl_el!=. & T4==1
					mat M[1,`z'] = e(b)
					mat M[1,`z'+1] = .

			 // regression coefficients
				regress atr_ch_bl_el T1 T3 T2 woredafixed*, vce(cluster kebele_code)
					mat P = r(table)
					mat N[1,`z'] = e(N)
					mat N[1,`z'+1] = .
					mat C1_beta[1,`z']=(_b[T1],_se[T1]\_b[T3],_se[T3]\_b[T2],_se[T2]\ ///
						.,.\.,.\.,.\.,.\.,.\.,.\.,.)
	
			// p-values
				forval x = 1/3 {
					mat C1_pvalues[`x',`z'] = P["pvalue",`x']
				}
				
			// test equality
				test T1 = T2
					mat C1_test1[1,`z'] = r(p),.
				test T2 = T3
					mat C1_test2[1,`z'] = r(p),.
				test T1 = T3
					mat C1_test3[1,`z'] = r(p),.

			// p-values
				local y = 1
				forval x = 11/13 {
					mat C1_pvalues[`x',`z'] = C1_test`y'[1,`z']
					local y = `y' + 1
				}
								   
			 // annotate arm coefficients with sig. stars
				forval x = 1/13 {
					mat C1_stars2[`x',`z'] = ((C1_pvalues[`x',`z']<=0.10) + (C1_pvalues[`x',`z']<=0.05) + (C1_pvalues[`x',`z']<=0.01)),0
				}
				
			mat drop P          // dropping matrix with regression results
	
	// now with baseline controls added 
	
	// run loop to create matrices as individual inputs for statmat
		local z = 3
		local x = 2 
			
			// means (matrix M)
				mean atr_ch_bl_el if atr_ch_bl_el!=. & T4==1
					mat M[1,`z'] = e(b)
					mat M[1,`z'+1] = .

			 // regression coefficients
				regress atr_ch_bl_el T1 T3 T2 `controls' m_bl_* woredafixed*, vce(cluster kebele_code)
					mat P = r(table)
					mat N[1,`z'] = e(N)
					mat N[1,`z'+1] = .
					mat C1_beta[1,`z']=(_b[T1],_se[T1]\_b[T3],_se[T3]\_b[T2],_se[T2]\ ///
						_b[hh_size_r1],_se[hh_size_r1]\_b[total_ch5_r1],_se[total_ch5_r1]\_b[log_tot_exp_ae],_se[log_tot_exp_ae]\ ///
						_b[hh_nearest_town_km],_se[hh_nearest_town_km]\_b[pf_marital_1_r1],_se[pf_marital_1_r1]\ ///
						_b[pf_age_r1],_se[pf_age_r1]\_b[pf_educ_r1],_se[pf_educ_r1])
	
			// p-values
				forval x = 1/10 {
					mat C1_pvalues[`x',`z'] = P["pvalue",`x']
				}
				
			// test equality
				test T1 = T2
					mat C1_test1[1,`z'] = r(p),.
				test T2 = T3
					mat C1_test2[1,`z'] = r(p),.
				test T1 = T3
					mat C1_test3[1,`z'] = r(p),.

			// p-values
				local y = 1
				forval x = 11/13 {
					mat C1_pvalues[`x',`z'] = C1_test`y'[1,`z']
					local y = `y' + 1
				}
								   
			 // annotate arm coefficients with sig. stars
				forval x = 1/13 {
					mat C1_stars2[`x',`z'] = ((C1_pvalues[`x',`z']<=0.10) + (C1_pvalues[`x',`z']<=0.05) + (C1_pvalues[`x',`z']<=0.01)),0
				}
				
			mat drop P          // dropping matrix with regression results
			
	// aggregate matrices into full statmat as input for frmttable         
		mat D2 = C1_beta\C1_test1\C1_test2\C1_test3\M\N	
		mat D = D1,D2
		mat C1_stars = C1_stars1,C1_stars2
		
	frmttable using "${tables}/all tables_${date}.doc", addtable landscape statmat(D) substat(1) annotate(C1_stars) asymbol(*,**,***) ///
		hlines(101000000000000000000010000010001) sdec(3\3\3\3\3\3\3\3\3\3\3\3\3\3\3\3\3\3\3\3\3\3\3\3\3\3\3\3\0\0) varlabels ///
		title("Annex Table A1: Attrition by treatment arm") nobl ///
		basefont(fs10) titlfont(fs13 b) ctitlfont(fs10) rtitlfont(fs10) notefont(fs9) colwidth(60 10 10 10 10 10 10 10 10 10 10) ///
		rtitle("EN with grants"\""\"EN without grants"\""\"CN with grants"\""\ ///
			"Household size"\""\"Number of children under the age of 5"\""\"Log of total monthly consumption expenditure per adult equivalent"\""\ ///
			"Distance to nearest town"\""\"Primary female is married monogamously"\""\"Age of primary female"\""\"Primary female has some education"\""\ ///
			"Test: EN with grants = CN with grants"\""\"Test: CN with grants = EN without grants"\""\"Test: EN with grants = EN without grants"\""\ ///
			"Mean of control"\""\"N") ///
		ctitles("","(1)","(2)","(3)","(4)"\"","Household is attrited from the endline sample","Household is attrited from the endline sample", ///
			"Household has no child in the endline anthropometry sample","Household has no child in the endline anthropometry sample") ///
		note("${footnote_attrit}") ///
		coljust(l;c)
	
	restore 
	
	//============================================================================//		
	}
		
	
	//==============================================================================//
	**# A2: BASELINE ANTHRO & SUBSEQUENT ATTRITION //
	{
	/*		use the approach from Fitzgerald Gotshalk and Moffit and run a regression of HAZ using baseline data with child age and gender, 
			mother's educ, a measure of HH assets, treatment group AND a dummy variable defined as 1 if the HH will be lost for any 
			reason before the endline. Unrest kebeles are excluded; woreda FE added.			*/
	
	preserve 
	drop if inlist(kebele,"Selamgae","Shimila albeta","Bad","Cherey")	// drop the 4 kebeles that weren't surveyed due to unrest (see the paper)
	keep if poultrycash_hh==1 
	
	// define matrices for statmat
		local i = 0
		local c : word count `rvars'
		local 3c = 6
		clear matrix
		mat	N				= J(1,`3c',0)		// sample N
		mat	M				= J(1,`3c',0)		// control means
		  
	// matrices for model with extended controls	
		mat C1_beta  		= J(8,`3c',0)		// coefficients for beta1 - beta3
		forval t = 1/3 {
			mat C1_test`t'	= J(1,`3c',0)		// tests of equality
		}
		mat C1_pvalues		= J(11,`3c',0)		// p-values
		mat C1_stars  		= J(11,`3c',0)		// annotated coefficients with ***
					 
		local controls hh_size_r1 total_ch5_r1 log_tot_exp_ae hh_nearest_town_km ///
			pf_marital_1_r1 pf_age_r1 pf_educ_r1
			
	// REGRESS HAZ ON HH LEVEL ATTRITION & T INDICATORS //
	
	// run loop to create matrices as individual inputs for statmat
		local z = 1
		local x = 1 
			
			// means (matrix M)
				mean ch_haz_r1 if T4==1
					mat M[1,`z'] = e(b)
					mat M[1,`z'+1] = .

			 // regression coefficients
				regress ch_haz_r1 T1 T3 T2 atr_hh_bl_el ch_agem_r1 ch_male_r1 pf_educ_r1 hh_totalland_asset_r1 woredafixed*, vce(cluster kebele_code)
					mat P = r(table)
					mat N[1,`z'] = e(N)
					mat N[1,`z'+1] = .
					mat C1_beta[1,`z']=(_b[T1],_se[T1]\_b[T3],_se[T3]\_b[T2],_se[T2]\ ///
						_b[atr_hh_bl_el],_se[atr_hh_bl_el]\_b[ch_agem_r1],_se[ch_agem_r1]\_b[ch_male_r1],_se[ch_male_r1]\ ///
						_b[pf_educ_r1],_se[pf_educ_r1]\_b[hh_totalland_asset_r1],_se[hh_totalland_asset_r1])
	
			// p-values
				forval x = 1/8 {
					mat C1_pvalues[`x',`z'] = P["pvalue",`x']
				}
				
			// test equality
				test T1 = T2
					mat C1_test1[1,`z'] = r(p),.
				test T2 = T3
					mat C1_test2[1,`z'] = r(p),.
				test T1 = T3
					mat C1_test3[1,`z'] = r(p),.

			// p-values
				local y = 1
				forval x = 9/11 {
					mat C1_pvalues[`x',`z'] = C1_test`y'[1,`z']
					local y = `y' + 1
				}
								   
			 // annotate arm coefficients with sig. stars
				forval x = 1/11 {
					mat C1_stars[`x',`z'] = ((C1_pvalues[`x',`z']<=0.10) + (C1_pvalues[`x',`z']<=0.05) + (C1_pvalues[`x',`z']<=0.01)),0
				}
				
			mat drop P          // dropping matrix with regression results
	
	// REGRESS HAZ ON CHILD LEVEL ATTRITION & T INDICATORS //
	
	// run loop to create matrices as individual inputs for statmat
		local z = 3
		local x = 2 
			
			// means (matrix M)
				mean ch_haz_r1 if T4==1
					mat M[1,`z'] = e(b)
					mat M[1,`z'+1] = .

			 // regression coefficients
				regress ch_haz_r1 T1 T3 T2 atr_ch_bl_el ch_agem_r1 ch_male_r1 pf_educ_r1 hh_totalland_asset_r1 woredafixed*
					mat P = r(table)
					mat N[1,`z'] = e(N)
					mat N[1,`z'+1] = .
					mat C1_beta[1,`z']=(_b[T1],_se[T1]\_b[T3],_se[T3]\_b[T2],_se[T2]\ ///
						_b[atr_ch_bl_el],_se[atr_ch_bl_el]\_b[ch_agem_r1],_se[ch_agem_r1]\_b[ch_male_r1],_se[ch_male_r1]\ ///
						_b[pf_educ_r1],_se[pf_educ_r1]\_b[hh_totalland_asset_r1],_se[hh_totalland_asset_r1])
	
			// p-values
				forval x = 1/8 {
					mat C1_pvalues[`x',`z'] = P["pvalue",`x']
				}
				
			// test equality
				test T1 = T2
					mat C1_test1[1,`z'] = r(p),.
				test T2 = T3
					mat C1_test2[1,`z'] = r(p),.
				test T1 = T3
					mat C1_test3[1,`z'] = r(p),.

			// p-values
				local y = 1
				forval x = 9/11 {
					mat C1_pvalues[`x',`z'] = C1_test`y'[1,`z']
					local y = `y' + 1
				}
								   
			 // annotate arm coefficients with sig. stars
				forval x = 1/11 {
					mat C1_stars[`x',`z'] = ((C1_pvalues[`x',`z']<=0.10) + (C1_pvalues[`x',`z']<=0.05) + (C1_pvalues[`x',`z']<=0.01)),0
				}
				
			mat drop P          // dropping matrix with regression results
			
	// REGRESS KNOWLEDGE ON HH LEVEL ATTRITION & T INDICATORS //
	
	// run loop to create matrices as individual inputs for statmat
		local z = 5
		local x = 3 
			
			// means (matrix M)
				mean knowledge_f_r1 if T4==1
					mat M[1,`z'] = e(b)
					mat M[1,`z'+1] = .

			 // regression coefficients
				regress knowledge_f_r1 T1 T3 T2 atr_hh_bl_el pf_educ_r1 hh_totalland_asset_r1 woredafixed*, vce(cluster kebele_code)
					mat P = r(table)
					mat N[1,`z'] = e(N)
					mat N[1,`z'+1] = .
					mat C1_beta[1,`z']=(_b[T1],_se[T1]\_b[T3],_se[T3]\_b[T2],_se[T2]\ ///
						_b[atr_hh_bl_el],_se[atr_hh_bl_el]\.,.\.,.\ ///
						_b[pf_educ_r1],_se[pf_educ_r1]\_b[hh_totalland_asset_r1],_se[hh_totalland_asset_r1])
	
			// p-values
				foreach x in 1 2 3 4 7 8 {
					mat C1_pvalues[`x',`z'] = P["pvalue",`x']
				}
				foreach x in 5 6 {
					mat C1_pvalues[`x',`z'] = .
				}
				
			// test equality
				test T1 = T2
					mat C1_test1[1,`z'] = r(p),.
				test T2 = T3
					mat C1_test2[1,`z'] = r(p),.
				test T1 = T3
					mat C1_test3[1,`z'] = r(p),.

			// p-values
				local y = 1
				forval x = 9/11 {
					mat C1_pvalues[`x',`z'] = C1_test`y'[1,`z']
					local y = `y' + 1
				}
								   
			 // annotate arm coefficients with sig. stars
				forval x = 1/11 {
					mat C1_stars[`x',`z'] = ((C1_pvalues[`x',`z']<=0.10) + (C1_pvalues[`x',`z']<=0.05) + (C1_pvalues[`x',`z']<=0.01)),0
				}
				
			mat drop P          // dropping matrix with regression results
			
	// aggregate matrices into full statmat as input for frmttable         
		mat D = C1_beta\C1_test1\C1_test2\C1_test3\M\N	
		
	frmttable using "${tables}/all tables_${date}.doc", addtable landscape statmat(D) substat(1) annotate(C1_stars) asymbol(*,**,***) ///
		hlines(10100000000000000010000010001) sdec(3\3\3\3\3\3\3\3\3\3\3\3\3\3\3\3\3\3\3\3\3\3\3\3\0\0) varlabels ///
		title("Annex Table A2: Baseline anthropometrics and subsequent attrition") nobl ///
		basefont(fs10) titlfont(fs13 b) ctitlfont(fs10) rtitlfont(fs10) notefont(fs9) colwidth(60 10 10 10 10 10 10 10 10 10 10) ///
		rtitle("EN with grants"\""\"EN without grants"\""\"CN with grants"\""\"Household attrited from the endline sample/"\"Household has no child in the endline anthropometry sample"\ ///
			"Child's age in months"\""\"Child is male"\""\"Mother has some education"\""\ ///
			"Household's asset and land index"\""\ ///
			"Test: EN with grants = CN with grants"\""\"Test: CN with grants = EN without grants"\""\"Test: EN with grants = EN without grants"\""\ ///
			"Mean of control"\""\"N") ///
		ctitles("","(1)","(2)","(3)"\"","Child's HAZ at baseline","Child's HAZ at baseline","Mother's IYCF knowledge score (0-7)") ///
		note("${footnote_attrit}") ///
		coljust(l;c)
		
	restore 
	}
	
	
	//==============================================================================//
	**# B1: SUMMARY BY OUTCOME FAMILIES //
	{
		
	// construct standardized variables 
	
		// mother's healthcare 
		preserve 
		tempfile temp_mo
			
			keep if !mi(j7_hew_hda_r3) & !mi(j7_hew_hda_home_r3) & !mi(j7_fdemo_3mo_r3) & ///
				!mi(j7_bcc_3mo_r3) & !mi(j7_cltsh_r3) & !mi(j7_odf_r3) & !mi(m2_topic_2a_f_r3)
			
			foreach var in j7_hew_hda j7_hew_hda_home j7_fdemo_3mo j7_bcc_3mo j7_cltsh j7_odf m2_topic_2a_f {
				
				// poor 
				qui sum `var'_r3 if T4==1 & poultrycash_hh==1
				local m = `r(mean)'
				local sd = `r(sd)'
				gen std_`var' = (`var'_r3 - `m')/`sd' if poultrycash_hh==1
			
				// less poor
				qui sum `var'_r3 if T4==1 & poultrycash_hh==0
				local m = `r(mean)'
				local sd = `r(sd)'
				replace std_`var' = (`var'_r3 - `m')/`sd' if poultrycash_hh==0
			}
			
			sum std_*
			egen std_healthcare_r3 = rowmean(std_*)	
			
			foreach var in j7_hew_hda j7_hew_hda_home j7_fdemo_3mo j7_bcc_3mo {
				
				// poor 
				qui sum `var'_r1 if T4==1 & poultrycash_hh==1
				local m = `r(mean)'
				local sd = `r(sd)'
				gen std_`var'_r1 = (`var'_r1 - `m')/`sd' if poultrycash_hh==1
			
				// less poor
				qui sum `var'_r1 if T4==1 & poultrycash_hh==0
				local m = `r(mean)'
				local sd = `r(sd)'
				replace std_`var'_r1 = (`var'_r1 - `m')/`sd' if poultrycash_hh==0
			}
			
			sum std_*r1
			egen std_healthcare_r1 = rowmean(std_*r1)	
			
			save `temp_mo', replace
			
		restore 
		
		// child's health 
		preserve 
		tempfile temp_ch
			
			keep if eic_age_r3 >= 6 & eic_age_r3 < 24
			keep if !mi(ch_vita_r3) & !mi(ch_weight_measure_r3) & !mi(ch_height_measure_r3) & !mi(ch_muac_measure_r3) & ///
				!mi(ch_underweight_r3) & !mi(ch_cpnp_r3)
			
			foreach var in ch_vita ch_weight_measure ch_height_measure ch_muac_measure ch_underweight ch_cpnp {
				
				// poor 
				qui sum `var'_r3 if T4==1 & poultrycash_hh==1
				local m = `r(mean)'
				local sd = `r(sd)'
				gen std_`var' = (`var'_r3 - `m')/`sd' if poultrycash_hh==1
			
				// less poor
				qui sum `var'_r3 if T4==1 & poultrycash_hh==0
				local m = `r(mean)'
				local sd = `r(sd)'
				replace std_`var' = (`var'_r3 - `m')/`sd' if poultrycash_hh==0
			}
			
			sum std_*
			egen std_childhealth_r3 = rowmean(std_*)
			save `temp_ch', replace
			
		restore 
		
		// merge in new indicators 
		merge 1:1 hh_id using `temp_mo', keepusing(std_healthcare_r3) nolabel
		drop _merge
		merge 1:1 hh_id using `temp_ch', keepusing(std_childhealth_r3) nolabel
		drop _merge			
		
	// TABLE 
		
	// define matrices for statmat
		local i = 0
		local c = 7
		local 2c = `c'*2
		clear matrix
		  
	// matrices for model with extended controls
		foreach sample in poor lpoor {
			mat	N_`sample'		= J(1,`2c',0)		// sample N
			mat	M_`sample'		= J(1,`2c',0)		// control means
			mat C1_beta_`sample' = J(2,`2c',0)		// coefficients for poor
		}
		forval t = 1/4 {
			mat C1_test`t'	= J(1,`2c',0)		// tests of equality
		}
		mat C1_pvalues		= J(12,`2c',0)		// p-values
		mat C1_stars  		= J(12,`2c',0)		// annotated coefficients with ***
		
		// model - anthro
			
			preserve
			
				use "$dataout/spir_child level data.dta", clear 
				
				keep if !mi(ch_haz_r3) & !mi(ch_stunted_r3) & !mi(ch_whz_r3) & !mi(ch_wasted_r3) & !mi(ch_waz_r3) & ///
					!mi(ch_underwght_r3) & !mi(ch_muac_r3) & !mi(ch_sam_r3)
				replace ch_haz_r3 = . if ch_agem_r3 < 6 | ch_agem_r3 >= 36
				replace ch_stunted_r3 = . if ch_agem_r3 < 6 | ch_agem_r3 >= 36
				
				local rvars ch_haz_r3 ch_stunted_r3
				local z = 1				
				forval i = 1/2 {    
					
					tokenize `rvars'
					local rnames "`rvars'"
					
					local end ``i''
					
					// EXTREMELY POOR SAMPLE
									
						// means
							mean `end' if poultrycash_hh==1 & T4==1
								mat M_poor[1,`z'] = e(b)
								mat M_poor[1,`z'+1] = .
								
						 // regression
							regress `end' T1_T3 T2 pf_educ_r1 ch_agem_r3 ch_male_r3 woredafixed* ///
								if poultrycash_hh==1, vce(cluster kebele)
								qui regress `end' T1_T3 T2 pf_educ_r1 ch_agem_r3 ch_male_r3 woredafixed* if poultrycash_hh==1
								estimates store poor
								mat P = r(table)
								mat N_poor[1,`z'] = e(N)
								mat N_poor[1,`z'+1] = .
								mat C1_beta_poor[1,`z']=(_b[T1_T3],_se[T1_T3]\_b[T2],_se[T2])
								
						// p-values
							forval x = 1/2 {
								mat C1_pvalues[`x',`z'] = P["pvalue",`x']
							}
							cap drop __0000*
						
						// test: EN vs CN
							test T1_T3 = T2
								mat C1_test1[1,`z'] = r(p),.
								
					// LESS POOR SAMPLE
							
						// means: LESS POOR
							mean `end' if poultrycash_hh==0 & T4==1
								mat M_lpoor[1,`z'] = e(b)
								mat M_lpoor[1,`z'+1] = .
								
						 // regression: LESS POOR
							regress `end' T1_T3 T2 pf_educ_r1 ch_agem_r3 ch_male_r3 woredafixed* ///
								if poultrycash_hh==0, vce(cluster kebele)
								qui regress `end' T1_T3 T2 pf_educ_r1 ch_agem_r3 ch_male_r3 woredafixed* if poultrycash_hh==0
								estimates store lpoor
								mat P = r(table)
								mat N_lpoor[1,`z'] = e(N)
								mat N_lpoor[1,`z'+1] = .
								mat C1_beta_lpoor[1,`z']=(_b[T1_T3],_se[T1_T3]\_b[T2],_se[T2])
				
						// p-values
							forval n = 1/2 {
								forval x = 5/6 {
									mat C1_pvalues[`x',`z'] = P["pvalue",`n']
								}
							}
						
						// test: EN = CN
							test T1_T3 = T2
								mat C1_test2[1,`z'] = r(p),.
								
					// TEST: poor vs less poor 
						qui suest poor lpoor, cluster(kebele)
							
						qui testnl [poor_mean]T1_T3 = [lpoor_mean]T1_T3
							mat C1_test3[1,`z'] = r(p),.
						qui testnl [poor_mean]T2 = [lpoor_mean]T2
							mat C1_test4[1,`z'] = r(p),.
						
					// p-values
						mat C1_pvalues[9,`z'] = C1_test1[1,`z']
						mat C1_pvalues[10,`z'] = C1_test2[1,`z']
						mat C1_pvalues[11,`z'] = C1_test3[1,`z']
						mat C1_pvalues[12,`z'] = C1_test4[1,`z']
						foreach n in 3 4 7 8 {
							mat C1_pvalues[`n',`z'] = .
						}
										   
					 // annotate arm coefficients with sig. stars 
						forval x = 1/12 {
							mat C1_stars[`x',`z'] = ((C1_pvalues[`x',`z']<=0.10) + (C1_pvalues[`x',`z']<=0.05) + (C1_pvalues[`x',`z']<=0.01)),0
						}
						
					mat drop P          // dropping matrix with regression results
					local z = `z'+2
					
				}     
			  
			restore
		
		// model - other outcomes 
			
			preserve 
			
				// restrict child's dietary diversity to the sample of 6-23mo olds 
					replace total_food_r3 = . if eic_age_r3 < 6 | eic_age_r3 >= 24
				
				local rvars std_healthcare std_childhealth knowledge_f wdds total_food
				
				local z = 5
				forval i = 1/5 {    
					
					tokenize `rvars'
					local end ``i''_r3
					capture confirm variable ``i''_r1
					if !_rc {
						local base 	``i''_r1  miss_``i''_r1
					}
					else {
						local base 	""
					}
					
					global controls pf_educ_r1 `base'
					if "`end'"=="std_childhealth_r3" {
						global controls pf_educ_r1 eic_age_r3 eic_male_r3
					}
					if "`end'"=="wdds_r3" {
						global controls pf_educ_r1 post_fast_r3 wdds_r1 miss_wdds_r1
					}
					if "`end'"=="total_food_r3" {
						global controls pf_educ_r1 eic_age_r3 eic_male_r3 post_fast_r3
					}
					
					di "OUTCOME: `end' ||| CONTROLS: $controls"
					
					// EXTREMELY POOR SAMPLE
									
						// means
							mean `end' if poultrycash_hh==1 & T4==1
								mat M_poor[1,`z'] = e(b)
								mat M_poor[1,`z'+1] = .
								
						 // regression
							regress `end' T1_T3 T2 $controls woredafixed* ///
								if poultrycash_hh==1, vce(cluster kebele_code)
								qui regress `end' T1_T3 T2 $controls woredafixed* if poultrycash_hh==1
								estimates store poor
								mat P = r(table)
								mat N_poor[1,`z'] = e(N)
								mat N_poor[1,`z'+1] = .
								mat C1_beta_poor[1,`z']=(_b[T1_T3],_se[T1_T3]\_b[T2],_se[T2])
								
						// p-values
							forval x = 1/2 {
								mat C1_pvalues[`x',`z'] = P["pvalue",`x']
							}
							cap drop __0000*
						
						// test: EN vs CN
							test T1_T3 = T2
								mat C1_test1[1,`z'] = r(p),.
								
					// LESS POOR SAMPLE
							
						// means: LESS POOR
							mean `end' if poultrycash_hh==0 & T4==1
								mat M_lpoor[1,`z'] = e(b)
								mat M_lpoor[1,`z'+1] = .
								
						 // regression: LESS POOR
							regress `end' T1_T3 T2 $controls woredafixed* ///
								if poultrycash_hh==0, vce(cluster kebele_code)
								qui regress `end' T1_T3 T2 $controls woredafixed* if poultrycash_hh==0
								estimates store lpoor
								mat P = r(table)
								mat N_lpoor[1,`z'] = e(N)
								mat N_lpoor[1,`z'+1] = .
								mat C1_beta_lpoor[1,`z']=(_b[T1_T3],_se[T1_T3]\_b[T2],_se[T2])
				
						// p-values
							forval n = 1/2 {
								forval x = 5/6 {
									mat C1_pvalues[`x',`z'] = P["pvalue",`n']
								}
							}
						
						// test: EN = CN
							test T1_T3 = T2
								mat C1_test2[1,`z'] = r(p),.
								
					// TEST: poor vs less poor 
						qui suest poor lpoor, cluster(kebele_code)
							
						qui testnl [poor_mean]T1_T3 = [lpoor_mean]T1_T3
							mat C1_test3[1,`z'] = r(p),.
						qui testnl [poor_mean]T2 = [lpoor_mean]T2
							mat C1_test4[1,`z'] = r(p),.
						
					// p-values
						mat C1_pvalues[9,`z'] = C1_test1[1,`z']
						mat C1_pvalues[10,`z'] = C1_test2[1,`z']
						mat C1_pvalues[11,`z'] = C1_test3[1,`z']
						mat C1_pvalues[12,`z'] = C1_test4[1,`z']
						foreach n in 3 4 7 8 {
							mat C1_pvalues[`n',`z'] = .
						}
										   
					 // annotate arm coefficients with sig. stars 
						forval x = 1/12 {
							mat C1_stars[`x',`z'] = ((C1_pvalues[`x',`z']<=0.10) + (C1_pvalues[`x',`z']<=0.05) + (C1_pvalues[`x',`z']<=0.01)),0
						}
						
					mat drop P          // dropping matrix with regression results
					local z = `z'+2
				  }     
				  
				restore 
		
		// aggregate matrices into full statmat as input for frmttable         
			mat I_0 = J(1,`c',.)
			mat D = C1_beta_poor\M_poor\N_poor\C1_beta_lpoor\M_lpoor\N_lpoor\C1_test1\C1_test2\C1_test3\C1_test4				
			mat list D
				
			frmttable using "${tables}/all tables_$date.doc", addtable landscape statmat(D) substat(1) annotate(C1_stars) asymbol(*,**,***) ///			
				hlines(101000000010000000100000001) sdec(3\3\3\3\3\3\0\0\3\3\3\3\3\3\0\0\3) ///
				title("Annex Table B1: Summary by family of outcome, extremely poor and less poor subsamples") nobl ///
				basefont(fs10) titlfont(fs13 b) ctitlfont(fs10) rtitlfont(fs10) notefont(fs9) colwidth(28 10 10 10 10 10 10 10 10) ///
				rtitle("Enhanced Nutrition (EN)"\""\"Core Nutrition (CN)"\""\"Mean of control"\""\"N"\""\ ///
					"Enhanced Nutrition (EN)"\""\"Core Nutrition (CN)"\""\"Mean of control"\""\"N"\""\ ///
					"Extremely poor: EN = CN"\""\"Less poor: EN = CN"\""\"EN: Extremely poor = less poor"\""\"CN: Extremely poor = less poor")  ///
				ctitles("","(1)","(2)","(3)","(4)","(5)","(6)","(7)"\"","Height-for-age z-score (HAZ)","Proportion stunted (HAZ<-2SD)", ///
					"ASTE: Mother's access to healthcare","ASTE: Child's health and services","Mother's IYCF knowledge score (0-7)", ///
					"Mother's Dietary Diversity Score (1-10)","Child's total food groups consumed (0-8)") ///
				note("Notes: Estimates from the DFSA SPIR endline survey sample. Columns (1) and (2) are restricted to children aged 6-35 months only, columns (4) and (7) are restricted to children aged 6-23 months only. Standard errors (in parentheses) are clustered at the kebele level. All models control for woreda level fixed effects and whether child's mother has any education; columns (1), (2), (4) and (7) additionally control for child's sex and age in months; columns (6) and (7) control for whether the day preceding the survey was a fast day. Asterisks indicate significance at the 10, 5 and 1 percent level and are calculated with respect to the standard errors.") ///
				coljust(l;c)

		cap mat drop D N M I_0 C1_beta C1_lincom1 C1_lincom2 C1_lincom3 C1_lincom4 C1_test1 C1_test2 C1_test3 C1_test4 C1_test5 C1_pvalues C1_stars
		
	}

	
	//==============================================================================//
	**# C1: ANTHROPOMETRY: 6-23mo //
	{
	preserve
	
		use "$datain/spir_child level data.dta", clear
	
	// run the model & output the table 
		keep if poultrycash_hh==1 
		keep if !mi(ch_haz_r3) & !mi(ch_stunted_r3) & !mi(ch_whz_r3) & !mi(ch_wasted_r3) & !mi(ch_waz_r3) & !mi(ch_underwght_r3) & !mi(ch_muac_r3) & !mi(ch_sam_r3)
		
	// define matrices for statmat
		local rvars ch_haz ch_stunted ch_whz ch_wasted ch_waz ch_underwght ch_muac ch_sam	
		local i = 0
		local c : word count `rvars'
		local 2c = `c'*2
		clear matrix
		mat	N				= J(1,`2c',0)		// sample N
		mat	M				= J(1,`2c',0)		// control means
		  
	// matrices for model with extended controls
		mat C1_beta  		= J(5,`2c',0)		// coefficients for beta1 - beta9
		forval t = 1/4 {
			mat C1_lincom`t'= J(1,`2c',0)		// effects for linear combinations
		}
		forval t = 1/4 {
			mat C1_test`t'	= J(1,`2c',0)		// tests of equality
		}
		mat C1_pvalues		= J(11,`2c',0)		// p-values
		mat C1_stars  		= J(11,`2c',0)		// annotated coefficients with ***
							 
	// run loop to create matrices as individual inputs for statmat
		local z = 1
		
		forval i = 1/`c' {    
			
			tokenize `rvars'
			local rnames "`rvars'"
			 
			local end ``i''_r3
			
			// means (matrix M)
				mean `end' if `end'!=. & T4==1 & ch_agem_r3 >= 6 & ch_agem_r3 < 24
					mat M[1,`z'] = e(b)
					mat M[1,`z'+1] = .

			 // regression coefficients
				/*	Poultry X EN				β1
					Cash X EN					β2
					Poultry X CN			 	β3
					Cash X CN					β4
					T3/EN w/o grants			β5		*/
				regress `end' $model_poor ch_agem_r3 ch_male_r3 pf_educ_r1 ///
						if ch_agem_r3 >= 6 & ch_agem_r3 < 24, vce(cluster kebele)
					mat P = r(table)
					mat N[1,`z'] = e(N)
					mat N[1,`z'+1] = .
					mat C1_beta[1,`z']=(_b[N_star_poultry],_se[N_star_poultry]\_b[N_star_cash],_se[N_star_cash]\ ///
						_b[T3],_se[T3]\_b[N_poultry],_se[N_poultry]\_b[N_cash],_se[N_cash])
			
			// p-values
				forval x = 1/5 {
					mat C1_pvalues[`x',`z'] = P["pvalue",`x']
				}
				
			// linear combination of betas
				/*	Linear combination of EN+grants		Weighted average: .5 * β1 + .5 * β2
					Linear combination of CN+grants		Weighted average: .5 * β4 + .5 * β5					*/					
				lincom N_star_poultry*$pr_wght_NS_P + N_star_cash*$pr_wght_NS_C
					mat C1_lincom1[1,`z'] = r(estimate), r(se)		// EN with grants
					mat C1_pvalues[6,`z'] = r(p)
				lincom N_poultry*$pr_wght_N_P + N_cash*$pr_wght_N_C
					mat C1_lincom2[1,`z'] = r(estimate), r(se)		// CN with grants
					mat C1_pvalues[7,`z'] = r(p)
	
			// test equality
				test N_star_poultry*$pr_wght_NS_P + N_star_cash*$pr_wght_NS_C = N_poultry*$pr_wght_N_P + N_cash*$pr_wght_N_C
					mat C1_test1[1,`z'] = r(p),.
				test N_star_poultry*$pr_wght_NS_P + N_star_cash*$pr_wght_NS_C = T3
					mat C1_test2[1,`z'] = r(p),.
				test N_star_poultry = N_star_cash
					mat C1_test3[1,`z'] = r(p),.
				test N_poultry = N_cash
					mat C1_test4[1,`z'] = r(p),.

			// p-values
				local y = 1
				forval x = 8/11 {
					mat C1_pvalues[`x',`z'] = C1_test`y'[1,`z']
					local y = `y' + 1
				}
								   
			 // annotate arm coefficients with sig. stars
			 // pvalues are stored in columns 1,2 of matrix C1_pvalues 
				forval x = 1/11 {
					mat C1_stars[`x',`z'] = ((C1_pvalues[`x',`z']<=0.10) + (C1_pvalues[`x',`z']<=0.05) + (C1_pvalues[`x',`z']<=0.01)),0
				}
			mat drop P          // dropping matrix with regression results
			local z = `z'+2
			
		}         	
	
	// aggregate matrices into full statmat as input for frmttable         
		mat D = C1_beta\C1_lincom1\C1_lincom2\C1_test1\C1_test2\C1_test3\C1_test4\M\N		
		mat list D
		
	frmttable using "${tables}/all tables_$date.doc", addtable landscape statmat(D) substat(1) annotate(C1_stars) asymbol(*,**,***) ///
		hlines($hlines) sdec($sdec) ///
		title("Annex Table C1: Anthropometrics, children 6-23 months of age") nobl ///
		basefont(fs10) titlfont(fs13 b) ctitlfont(fs10) rtitlfont(fs10) notefont(fs9) colwidth(28 10 10 10 10 10 10 10 10 10 10) ///
		rtitle("EN x Poultry"\""\"EN x Cash"\""\"EN without grants"\""\"CN x Poultry"\""\"CN x Cash"\""\ ///
			"Average effect of EN with grants"\""\"Average effect of CN with grants"\""\ ///
			"Test: EN with grants = CN with grants"\""\"Test: EN with grants = EN without grants"\""\"Test: EN x Poultry = EN x Cash"\""\"Test: CN x Poultry = CN x Cash"\""\"Mean of control"\""\"N")  ///
		ctitles("","(1)","(2)","(3)","(4)","(5)","(6)","(7)","(8)"\"","Height-for-age z-score (HAZ)","Proportion stunted (HAZ<-2SD)","Weight-for-height z-score (WHZ)", ///
			"Proportion wasted (WHZ<-2SD)","Weight-for-age z-score (WAZ)","Proportion underweight (WAZ<-2SD)","Mid-upper arm circumference (MUAC)", ///
			"Proportion of severe acute malnutrition (MUAC<11.5 cm)") ///
		note("Notes: Estimates from the DFSA SPIR endline survey for children of age 6 to 35 months from the subsample of extremely poor households. Standard errors (in parentheses) are clustered at the kebele level. All models control for woreda level fixed effects, child's sex and age in months, and whether child's mother has any education. Asterisks indicate significance at the 10, 5 and 1 percent level and are calculated with respect to the standard errors. ${footnote_poor}") ///
		coljust(l;c)	
			
	cap mat drop D N M I_0 C1_beta C1_lincom1 C1_lincom2 C1_lincom3 C1_lincom4 C1_test1 C1_test2 C1_test3 C1_test4 C1_test5 C1_test6 C1_test7 C1_pvalues C1_stars
	restore
	}
	

	//==============================================================================//
	**# C2: ANTHROPOMETRY: 24-35mo //
	{
	preserve
	
		use "$datain/spir_child level data.dta", clear
	
	// run the model & output the table 
		keep if poultrycash_hh==1 
		keep if !mi(ch_haz_r3) & !mi(ch_stunted_r3) & !mi(ch_whz_r3) & !mi(ch_wasted_r3) & !mi(ch_waz_r3) & !mi(ch_underwght_r3) & !mi(ch_muac_r3) & !mi(ch_sam_r3)
		
	// define matrices for statmat
		local rvars ch_haz ch_stunted ch_whz ch_wasted ch_waz ch_underwght ch_muac ch_sam	
		local i = 0
		local c : word count `rvars'
		local 2c = `c'*2
		clear matrix
		mat	N				= J(1,`2c',0)		// sample N
		mat	M				= J(1,`2c',0)		// control means
		  
	// matrices for model with extended controls
		mat C1_beta  		= J(5,`2c',0)		// coefficients for beta1 - beta9
		forval t = 1/4 {
			mat C1_lincom`t'= J(1,`2c',0)		// effects for linear combinations
		}
		forval t = 1/4 {
			mat C1_test`t'	= J(1,`2c',0)		// tests of equality
		}
		mat C1_pvalues		= J(11,`2c',0)		// p-values
		mat C1_stars  		= J(11,`2c',0)		// annotated coefficients with ***
							 
	// run loop to create matrices as individual inputs for statmat
		local z = 1
		
		forval i = 1/`c' {    
			
			tokenize `rvars'
			local rnames "`rvars'"
			 
			local end ``i''_r3
			
			// means (matrix M)
				mean `end' if `end'!=. & T4==1 & ch_agem_r3 >= 24 & ch_agem_r3 < 36
					mat M[1,`z'] = e(b)
					mat M[1,`z'+1] = .

			 // regression coefficients
				/*	Poultry X EN				β1
					Cash X EN					β2
					Poultry X CN			 	β3
					Cash X CN					β4
					T3/EN w/o grants			β5		*/
				regress `end' $model_poor ch_agem_r3 ch_male_r3 pf_educ_r1 ///
						if ch_agem_r3 >= 24 & ch_agem_r3 < 36, vce(cluster kebele)
					mat P = r(table)
					mat N[1,`z'] = e(N)
					mat N[1,`z'+1] = .
					mat C1_beta[1,`z']=(_b[N_star_poultry],_se[N_star_poultry]\_b[N_star_cash],_se[N_star_cash]\ ///
						_b[T3],_se[T3]\_b[N_poultry],_se[N_poultry]\_b[N_cash],_se[N_cash])
			
			// p-values
				forval x = 1/5 {
					mat C1_pvalues[`x',`z'] = P["pvalue",`x']
				}
				
			// linear combination of betas
				/*	Linear combination of EN+grants		Weighted average: .5 * β1 + .5 * β2
					Linear combination of CN+grants		Weighted average: .5 * β4 + .5 * β5					*/					
				lincom N_star_poultry*$pr_wght_NS_P + N_star_cash*$pr_wght_NS_C
					mat C1_lincom1[1,`z'] = r(estimate), r(se)		// EN with grants
					mat C1_pvalues[6,`z'] = r(p)
				lincom N_poultry*$pr_wght_N_P + N_cash*$pr_wght_N_C
					mat C1_lincom2[1,`z'] = r(estimate), r(se)		// CN with grants
					mat C1_pvalues[7,`z'] = r(p)
	
			// test equality
				test N_star_poultry*$pr_wght_NS_P + N_star_cash*$pr_wght_NS_C = N_poultry*$pr_wght_N_P + N_cash*$pr_wght_N_C
					mat C1_test1[1,`z'] = r(p),.
				test N_star_poultry*$pr_wght_NS_P + N_star_cash*$pr_wght_NS_C = T3
					mat C1_test2[1,`z'] = r(p),.
				test N_star_poultry = N_star_cash
					mat C1_test3[1,`z'] = r(p),.
				test N_poultry = N_cash
					mat C1_test4[1,`z'] = r(p),.

			// p-values
				local y = 1
				forval x = 8/11 {
					mat C1_pvalues[`x',`z'] = C1_test`y'[1,`z']
					local y = `y' + 1
				}
								   
			 // annotate arm coefficients with sig. stars
			 // pvalues are stored in columns 1,2 of matrix C1_pvalues 
				forval x = 1/11 {
					mat C1_stars[`x',`z'] = ((C1_pvalues[`x',`z']<=0.10) + (C1_pvalues[`x',`z']<=0.05) + (C1_pvalues[`x',`z']<=0.01)),0
				}
			mat drop P          // dropping matrix with regression results
			local z = `z'+2
			
		}         	
	
	// aggregate matrices into full statmat as input for frmttable         
		mat D = C1_beta\C1_lincom1\C1_lincom2\C1_test1\C1_test2\C1_test3\C1_test4\M\N		
		mat list D
		
	frmttable using "${tables}/all tables_$date.doc", addtable landscape statmat(D) substat(1) annotate(C1_stars) asymbol(*,**,***) ///
		hlines($hlines) sdec($sdec) ///
		title("Annex Table C2: Anthropometrics, children 24-35 months of age") nobl ///
		basefont(fs10) titlfont(fs13 b) ctitlfont(fs10) rtitlfont(fs10) notefont(fs9) colwidth(28 10 10 10 10 10 10 10 10 10 10) ///
		rtitle("EN x Poultry"\""\"EN x Cash"\""\"EN without grants"\""\"CN x Poultry"\""\"CN x Cash"\""\ ///
			"Average effect of EN with grants"\""\"Average effect of CN with grants"\""\ ///
			"Test: EN with grants = CN with grants"\""\"Test: EN with grants = EN without grants"\""\"Test: EN x Poultry = EN x Cash"\""\"Test: CN x Poultry = CN x Cash"\""\"Mean of control"\""\"N")  ///
		ctitles("","(1)","(2)","(3)","(4)","(5)","(6)","(7)","(8)"\"","Height-for-age z-score (HAZ)","Proportion stunted (HAZ<-2SD)","Weight-for-height z-score (WHZ)", ///
			"Proportion wasted (WHZ<-2SD)","Weight-for-age z-score (WAZ)","Proportion underweight (WAZ<-2SD)","Mid-upper arm circumference (MUAC)", ///
			"Proportion of severe acute malnutrition (MUAC<11.5 cm)") ///
		note("Notes: Estimates from the DFSA SPIR endline survey for children of age 6 to 35 months from the subsample of extremely poor households. Standard errors (in parentheses) are clustered at the kebele level. All models control for woreda level fixed effects, child's sex and age in months, and whether child's mother has any education. Asterisks indicate significance at the 10, 5 and 1 percent level and are calculated with respect to the standard errors. ${footnote_poor}") ///
		coljust(l;c)	
			
	cap mat drop D N M I_0 C1_beta C1_lincom1 C1_lincom2 C1_lincom3 C1_lincom4 C1_test1 C1_test2 C1_test3 C1_test4 C1_test5 C1_test6 C1_test7 C1_pvalues C1_stars
	restore
	}
	
	
	//==============================================================================//
	**# D1: PERSISTENCE OF FOOD POLICY 2022 RESULTS //
	{
	
	// MIDLINE 
	
	preserve
	
	use "$datain/spir_midline data.dta", clear 
	
	// labels
		local lab1 : variable label foodgroup5_r2
		local lab2 : variable label wfoodgroup6_r2
		local lab3 : variable label poultry_owned_m_r2
	
	// means
		sum foodgroup5_r2 if (ch_agem_r2 > 5 & ch_agem_r2<=23) & T4==1
		local m1 = round(`r(mean)',0.001)
		sum wfoodgroup6_r2 if T4==1
		local m2 = round(`r(mean)',0.001)
		sum poultry_owned_m_r2 if T4==1
		local m3 = round(`r(mean)',0.001)

		outreg, clear
	
	// TABLE 2 COL 5 [MIDLINE]
		reg foodgroup5_r2 poultry_dummy cash_dummy L_star_other N_star N_star_poultry N_star_cash pred_consumpt_bl post_fast_r2 ///
			pf_educ_r1 pm_educ_r1 hh_size_r1 ///
			ch_agem_r2 ch_male_r2 mean_ch miss_pf_educ_r1 miss_pm_educ_r1 woredafixed* if (ch_agem_r2 > 5 & ch_agem_r2<=23), ///
			vce(cluster kebele_code)
		outreg, merge bdec(3) se starloc(1) starlevels(10 5 1) varlabels ctitle("","(1)" \ "","`lab1'") ///
			drop(miss_pf_educ_r1 miss_pm_educ_r1 woredafixed*) 
			
	// TABLE 2 COL 6 [MIDLINE]
		reg wfoodgroup6_r2 poultry_dummy cash_dummy L_star_other N_star N_star_poultry N_star_cash pred_consumpt_bl post_fast_r2 ///
			pf_educ_r1 pm_educ_r1 hh_size_r1 ///
			mean_w miss_pf_educ_r1 miss_pm_educ_r1 woredafixed* if n==1, vce(cluster kebele_code)
		outreg, merge bdec(3) se starloc(1) starlevels(10 5 1) varlabels ctitle("","(2)" \ "","`lab2'") ///
			drop(miss_pf_educ_r1 miss_pm_educ_r1 woredafixed*) 
		
	// TABLE 4 COL 2 [MIDLINE]
		reg poultry_owned_m_r2 poultry_dummy cash_dummy L_star_other N_star pred_consumpt_bl pf_educ_r1 pm_educ_r1 hh_size_r1 miss_pf_educ_r1 miss_pm_educ_r1 ///
			poultry_owned_m_r1 woredafixed* if n==1, vce(cluster kebele_code)
		outreg, merge bdec(3) se starloc(1) starlevels(10 5 1) varlabels ctitle("","(3)" \ "","`lab3'") ///
			drop(miss_pf_educ_r1 miss_pm_educ_r1 woredafixed*) 
		
	// TABLE 5 COL 1 [MIDLINE]
		reg foodgroup5_r2 poultry_owned_m_r2 N_star pred_consumpt_bl post_fast_r2 pf_educ_r1 pm_educ_r1 hh_size_r1 ch_agem_r2 ch_male_r2 ///
			mean_ch miss_pf_educ_r1 miss_pm_educ_r1 woredafixed* if ch_agem_r2>5 & ch_agem_r2<=23, vce(cluster kebele_code)
		outreg, merge bdec(3) se starloc(1) starlevels(10 5 1) varlabels ctitle("","(4)" \ "","`lab1'") ///
			drop(miss_pf_educ_r1 miss_pm_educ_r1 woredafixed*) 
	
	// TABLE 5 COL 2 [MIDLINE]
		reg wfoodgroup6_r2 poultry_owned_m_r2 N_star pred_consumpt_bl post_fast_r2 pf_educ_r1 pm_educ_r1 hh_size_r1 ///
			mean_w miss_pf_educ_r1 miss_pm_educ_r1 woredafixed* if n==1, vce(cluster kebele_code)
		outreg, merge bdec(3) se starloc(1) starlevels(10 5 1) varlabels ctitle("","(5)" \ "","`lab2'") ///
			drop(miss_pf_educ_r1 miss_pm_educ_r1 woredafixed*)
	
	restore 
	
	// ENDLINE 
	
	preserve
	
	// replace missing BL outcomes with 0s + control for it
		foreach var in pf_educ pm_educ {
			gen miss_`var'_r1 = missing(`var'_r1)
			replace `var'_r1 = 0 if missing(`var'_r1)
		}	
				
	// means
		qui sum foodgroup5_r3 if (eic_age_r3 > 5 & eic_age_r3<=23) & T4==1
		local m6 = round(`r(mean)',0.001)
		qui sum wfoodgroup6_r3 if T4==1
		local m7 = round(`r(mean)',0.001)
		qui sum poultry_owned_m_r3 if T4==1
		local m8 = round(`r(mean)',0.001)
				
	// TABLE 2 COL 5 [ENDLINE]
		reg foodgroup5_r3 poultry_dummy cash_dummy L_star_other T1_T3 N_star_poultry N_star_cash post_fast_r3 pf_educ_r1 pm_educ_r1 hh_size_r1 ///
			eic_age_r3 eic_male_r3 miss_pf_educ_r1 miss_pm_educ_r1 woredafixed* if (eic_age_r3 > 5 & eic_age_r3<=23), vce(cluster kebele_code)
		outreg, merge bdec(3) se starloc(1) starlevels(10 5 1) varlabels ctitle("","(6)" \ "","`lab1'") ///
			drop(miss_pf_educ_r1 miss_pm_educ_r1 woredafixed*) 
			
	// TABLE 2 COL 6 [ENDLINE]
		reg wfoodgroup6_r3 poultry_dummy cash_dummy L_star_other T1_T3 N_star_poultry N_star_cash post_fast_r3 pf_educ_r1 pm_educ_r1 hh_size_r1 ///
			miss_pf_educ_r1 miss_pm_educ_r1 woredafixed*, vce(cluster kebele_code)
		outreg, merge bdec(3) se starloc(1) starlevels(10 5 1) varlabels ctitle("","(7)" \ "","`lab2'") ///
			drop(miss_pf_educ_r1 miss_pm_educ_r1 woredafixed*) 
		
	// TABLE 4 COL 2 [ENDLINE]
		reg poultry_owned_m_r3 poultry_dummy cash_dummy L_star_other T1_T3 pf_educ_r1 pm_educ_r1 hh_size_r1 miss_pf_educ_r1 miss_pm_educ_r1 ///
			poultry_owned_m_r1 woredafixed*, vce(cluster kebele_code)
		outreg, merge bdec(3) se starloc(1) starlevels(10 5 1) varlabels ctitle("","(8)" \ "","`lab3'") ///
			drop(miss_pf_educ_r1 miss_pm_educ_r1 woredafixed*) 
		
	// TABLE 5 COL 1 [ENDLINE]
		reg foodgroup5_r3 poultry_owned_m_r3 T1_T3 post_fast_r3 pf_educ_r1 pm_educ_r1 hh_size_r1 eic_age_r3 eic_male_r3 ///
			miss_pf_educ_r1 miss_pm_educ_r1 woredafixed* if eic_age_r3>5 & eic_age_r3<=23, vce(cluster kebele_code)
		outreg, merge bdec(3) se starloc(1) starlevels(10 5 1) varlabels ctitle("","(9)" \ "","`lab1'") ///
			drop(miss_pf_educ_r1 miss_pm_educ_r1 woredafixed*) 
	
	// TABLE 5 COL 2 [ENDLINE]
		reg wfoodgroup6_r3 poultry_owned_m_r3 T1_T3 post_fast_r3 pf_educ_r1 pm_educ_r1 hh_size_r1 ///
			miss_pf_educ_r1 miss_pm_educ_r1 woredafixed*, vce(cluster kebele_code)
		outreg, merge bdec(3) se starloc(1) starlevels(10 5 1) varlabels ctitle("","(10)" \ "","`lab2'") ///
			drop(miss_pf_educ_r1 miss_pm_educ_r1 woredafixed*)
		
	// table for full sample	
		outreg using "${tables}/all tables_$date.doc", landscape varlabels replay addtable ///
			title("Annex Table D1: Testing the persistence of results of poultry ownership: midline and endline") ///
			coljust(l;c) colwidth(60 8 8 8 8 8 8 8 8 8 8) basefont(fs10) titlfont(fs13 b) ///
			addrows("Mean of control","`m1'","`m2'","`m3'","`m1'","`m2'","`m6'","`m7'","`m8'","`m6'","`m7'") ///
			note("Notes: Estimates from the DFSA SPIR endline and midline survey samples; the midline results are from Alderman et al 2022, Food Policy. Standard errors (in parentheses) are clustered at the kebele level. All models control for woreda level fixed effects. Asterisks indicate significance at the 10, 5 and 1 percent level and are calculated with respect to the standard errors. Where primary female's/male's education is missing, it is replaced with 0 and respectively controlled for.") 
	
	restore
	
	}

}
